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Notation 


A  -  Plant  matrix 
B  -  Control  matrix 
C  -  Output  matrix 

D  -  Nodal  attitude-evaluated  actuator  location  matrix 
E  -  Damping  matrix 

6  -  Optimal  steady  full -state  feedback  gain  matrix 
g  -  Generalized  coordinate 
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T  -  Quadratic  performance  index 
K  -  Suboptimal  feedback  gain  matrix;  stiffness  matrix 

K  -  Untransformed  feedback  gain  matrix 

K  -  Design  feedback  gain  matrix 

M  -  Mass  matrix 

n,  -  Number  of  actuators 

a 

n  -  Number  of  modes 
m 

Q  -  State  weighting  matrix 
R  -  Control  weighting  matrix 

S  -  Singular  value  matrix;  solution  to  steady  state  matrix  Riccati 
equation 

T  -  Transformation  matrix 
u  -  Control  vector 
V  -  Matrix  of  right  singular  vectors 
W  -  Matrix  of  left  singular  vectors 
w  -  Transformed  output  vector 
x  -  State  vector 
y  -  Output  vector 

v  -  Transformed  control  vector 
r  -  Transformation  matrix 

<  -  Sum  of  K.j  's 

A  -  Sum  of  K.. '  s 
4>  -  Modal  matrix 

¥  -  Structural  mode  amplitude  matrix 
ip  -  Element  of  ¥  matrix 
a  -  Singular  value 

E  -  Suranatlon  when  used  with  indices;  singular  value  matrix 
C  -  Damping  ratio 
(u  -  Natural  frequency 
Q  -  w2 
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Notation  (Continued) 


Subscripts 


Critical  or  controlled  modes 
Index  for  variable 
Index  for  variable 
Index  for  variable 

Partition  of  V  matrix  associated  with  zero  singular  values; 
position  elements 

Partition  of  W  matrix  associated  with  zero  singular  values; 
residual  modes 

Partition  of  V  or  W  matrices  associated  with  non-zero  singular 
values;  suppressed  modes 
Unmodeled  modes 
Velocity  elements 
N-l  Control  matrices 


Superscripts 


*  -  Transformed  matrix 

+  -  Inverse  of  the  singular  value  decomposition  of  a  matrix 
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Abstract 

Direct  output  feedback  control  methods  are  used  to  develop  a  multiple- 
input  multiple-output  controller.  The  controller  is  then  applied  to  the 
Charles  Stark  Draper  Laboratory  2  (CSDL  2)  model.  The  CSDL  2  model  is  a 
sophisticated  optical  space  structure  representative  of  large  flexible 
space  structures.  This  model  consists  of  59  nodes  and  23  lumped  masses. 

The  beam  elements  are  fully  connected  and  may  support  axial,  transverse, 
and  torsional  deformations.  NASTRAN  is  employed  to  generate  modal  approx¬ 
imations  of  the  model  as  well  as  the  mode  shapes  and  frequencies  of  the 
resulting  modes.  Of  the  numerous  modes  available  for  the  model,  only  the 
first  36  modes  are  utilized  and  implemented  in  the  controller. 

The  control  problem  is  formulated  in  state  space  form  and  direct 
output  feedback  is  implemented.  The  state  is  represented  as  modal  ampli¬ 
tudes  and  rates.  System  outputs  are  obtained  by  rate  sensors  and  control 
is  applied  by  point  force  actua tors. vSi nee  the  observation  matrices  are 
not  of  full  rank,  their  generalized  (Penrose)  inverses  are  obtained.  These 
and  the  steady  full -state  feedback  optimal  control  gains  are  used  to  de¬ 
termine  suboptimal  feedback  gains.  In  the  multiple  controller,  individual 
controllers  actively  control  the  modes  assigned  to  it,  but  the  modes  as¬ 
signed  to  the  other  controllers  interact  with  all  controllers  causing  con¬ 
trol  and  observation  spillover.  The  spillover  that  appears  in  coupled  de¬ 
centralized  systems  is  eliminated  using  transformation  matrices  obtained 
via  modal  suppression  techniques.  The  transformation  matrices  are  applied 
to  the  feedback  gains  to  drive  the  spillover  terms  to  zero  and  decouple 
the  controllers. 


viii 


Control  was  accomplished  using  four  controllers  without  residual  modes 
and  three  controllers  with  residual  modes.  Conditions  for  which  the  sta¬ 
bility  of  the  model  is  assured  are  developed.  Arrangement  of  modal  assign 
ment  to  controllers  was  varied  to  examine  modal  selection  and  grouping. 
Sensitivity  to  uncontrolled  residual  modes  was  examined.  Direct  output 
feedback  response  was  compared  to  that  obtained  using  optimal  full -state 
feedback.  Full  controller  decoupling  was  achieved  and  stability  of  the 
non-rigid  body  modes  maintained.  The  investigation  demonstrated  the  fea¬ 
sibility  of  using  decentralized  direct  output  feedback  in  maintaining  sys¬ 
tem  stability. 


DECENTRALIZED  CONTROL  OF  A  LARGE  SPACE 
STRUCTURE  USING  DIRECT  OUTPUT  FEEDBACK 

I.  Introduction 

The  study  of  Large  Space  Structure  (LSS)  Control  is  a  growing  field 
of  interest  primarily  due  to  the  involvement  by  the  Department  of  Defense 
(DoD)  and  the  National  Aeronautics  and  Space  Administration  (15,16).  Tech 
nology  and  politics  are  leading  us  into  space  to  perform  scientific,  ex¬ 
ploratory,  and  eventually  operational  missions  that  use  large  structures. 
The  future  of  space  is  going  to  see  deployment  of  large  space  structures 
for  a  wide  spectrum  of  missions  including  communications,  space  defense, 
power  generation,  manufacturing,  and  research.  Future  DoD  large  space  sys 
tern  concepts  such  as  High  Altitude  Large  Optics  (HALO),  High  Energy  Laser 
Optics  (HELO),  and  Millimeter  (MM)  may  be  among  those  deployed.  These  sys 
terns  combine  large  size  with  stringent  line-of-sight  (LOS)  pointing  and 
jitter  requirements.  Coupled  with  these  performance  requirements  are  the 
needs  to  minimize  spacecraft  mass.  In  general,  the  constraints  on  large 
space  structures  result  in  highly  flexible  spacecraft  for  which  passive 
methods  of  material  selection  and  construction  techniques  do  not  provide 
adequate  control.  This  observation  has  led  to  the  idea  of  building  con¬ 
trol-configured  spacecraft  using  active  stability-augmentation  systems  in¬ 
tegrated  with  the  space  structure. 

To  achieve  this  control -configured  spacecraft,  the  control  designer 
must  understand  the  LSS  control  problem.  Some  characteristics  of  the  LSS 
control  problem  are:  (1)  infinite  dimensional  in  theory,  i.e.  distributed 


parameter  system,  and  very  large  dimensional  in  practice,  (2)  they  have 
many  low  resonant  frequencies  and  often  these  appear  in  closely-spaced 
"clumps",  (3)  their  natural  damping  is  very  light,  (4)  prediction  of  their 
behavior  in  space  by  ground  testing  is  quite  limited,  and  (5)  requirements 
for  shape,  orientation,  alignment,  vibration  suppression,  and  pointing  ac¬ 
curacy  are  extremely  rigorous  (6,18).  The  first  characteristic  is  consid¬ 
ered  fundamental  and  many  control  engineers  are  developing  and  using  ratio¬ 
nal  methods  to  reduce  theoretical  models  to  finite  dimensional  ones  that 
are  representative  of  the  actual  system  and  enable  control  systems  to  de¬ 
liver  adequate  performance  (4,6,9,15,19,21).  The  LSS  is  modeled  as  a  re¬ 
duced  order  discrete  system  via  finite  element  modeling  that  results  in  a 
finite-dimensional  modal  representation.  The  number  of  modeled  modes  in 
the  physical  model  is  referred  to  as  the  order  of  the  system.  When  a  dis¬ 
tributed  system  is  "discretized",  the  resulting  physical  model  suffers  er¬ 
rors  that  can  result  in  overall  system  instability  when  not  accounted  for 
(1,10).  From  this  physical  model  a  reduced  order  evaluation  model  is  se¬ 
lected  to  maintain  model  fidelity.  Since  the  evaluation  model  may  contain 
more  modes  than  can  be  controlled,  a  design  model  is  determined  to  include 
modes  which  if  ignored  would  degrade  system  performance  beyond  mission  re¬ 
quirements  (18).  Critical  modes,  denoted  x  ,  are  those  modes  of  the  design 
model  which  are  chosen  to  be  explicitly  controlled  in  order  to  assure  sta¬ 
bility  and  achieve  performance  requirements  for  the  system.  Residual  modes, 
denoted  xr,  are  those  modes  which  exist  in  the  infinite  dimensional  system 
that  are  not  critical  in  the  sense  defined  above.  Observation  spillover  is 
the  contamination  of  the  sensor  outputs  by  modes  not  contained  in  a  control¬ 
ler.  Control  spillover  is  the  excitation  of  mode  dynamics  due  to  a  control¬ 
ler  that  does  not  contain  those  modes.  Residual  modes  may  be  subdivided 


into  design,  evaluation,  physical,  and  unmodeled  modes.  Design  residual 

modes  are  those  which  are  included  in  the  controller  design  model,  but  not 

explicitly  controlled.  Evaluation  residual  modes  are  those  modes  included 

in  the  evaluation  model  but  not  in  the  design  model.  Physical  residual  modes 

are  those  modes  which  are  included  in  the  physical  model,  but  not  in  the 

evaluation  model.  Unmodeled  residual  modes,  denoted  x  .  are  those  modes 

urn 

which  exist,  but  are  not  included  in  any  finite-dimensional  model  of  the 
system  under  study  (18). 

The  specific  mission  of  the  LSS  forces  the  performance  requirements 
which  can  only  be  met  through  active  control.  There  are  many  appropriate 
control  strategies  for  large  flexible  space  structures  (1,2,10,18,20).  The 
choices  center  on  the  theme  of  constant-gain  linear  full  state  feedback 
control.  Of  the  choices  available,  decentralized  control  offers  several 
advantages.  In  decentralized  control  each  controller  is  designed  to  con¬ 
trol  some  of  the  critical  modes  allowing  the  use  of  low-order  computer  al¬ 
gorithms.  The  feedback  system  structure  used,  shown  in  Figure  1,  provides 
each  control  channel  with  the  same  information  structure.  This  method  of 
decentralized  control  using  modern  (cost  optimal)  control  for  finite  ele¬ 
ment  models  of  LSS  has  been  demonstrated  by  several  researchers  (1,2,13). 

They  found  transformation  matrices  applied  to  feedback  gains  proved  an  ef¬ 
fective  method  for  eliminating  spillover  and  enhancing  system  performance. 
Optimal  control  methods  which  require  estimators  to  reconstruct  the  states, 
however,  require  four  states  for  each  critical  mode  (mode  amplitude  and 
rate  and  estimated  amplitude  and  rate).  With  four  states  per  mode  the  com¬ 
putation  burden  rises  rapidly  with  the  number  of  critical  modes,  especially 
computation  of  the  algebraic  Riccati  equations  (8,11,20).  This  has  moti¬ 
vated  the  use  of  suboptimal  techniques  that  do  not  require  an  observer 


Figure  1.  Decentralized  Output  Feedback  Structure 

and  consequently  either  greatly  reduces  the  computational  burden  for  con¬ 
trolling  the  same  number  of  modes  or  allows  control  of  twice  as  many  crit¬ 
ical  modes.  This  concept  of  a  decentralized  direct  output  feedback  (DOFB) 
controller  was  first  proposed  by  Calico  (6)  and  will  be  discussed  in  de¬ 
tail  in  Section  III. 

Even  with  an  expanded  number  of  controllers  made  available  by  decen¬ 
tralized  control,  the  number  of  modes  that  may  be  controlled  is  small  com¬ 
pared  to  the  number  of  modes  that  exist  for  a  given  structure.  Consequently 
the  selection  of  modes  to  be  controlled  must  be  made  carefully.  Only  those 
modes  affecting  performance  need  be  controlled.  The  terms  "controlled" 
and  "critical"  will  be  used  interchangeably  for  these  modes.  The  remaining 
uncontrolled  modes  fall  into  three  categories:  suppressed,  residual  and 


unmodeled.  Natural  damping  in  the  structure  will  normally  prevent  insta¬ 
bilities  arising  from  the  higher  frequency  modes,  so,  for  model  simplicity, 
those  are  truncated  and  left  as  unmodeled  modes.  The  remaining  uncontrol¬ 
led  modes  are  modeled  modes.  Of  these,  some  may  have  destabilizing  effects 
due  to  spillover  and  therefore  have  to  be  made  transparent  to  the  control¬ 
ler.  These  are  referred  to  as  suppressed  modes.  The  last  mode  group  is 
modeled,  uncontrolled  and  unsuppressed.  These  are  the  residual  modes  and 
may  move  freely  when  control  is  applied.  They  may  become  more  stable,  less 
stable  or  unstable  due  to  control  and  observation  spillover  from  the  crit¬ 
ical  modes. 

The  intent  of  this  thesis  is  to  apply  direct  output  feedback  in  devel¬ 
oping  a  decentralized  control  system  of  three  or  four  controllers.  This 
control  system  will  be  applied  to  the  Charles  Stark  Draper  Laboratory,  Inc. 
model  #2  (CSDL  2).  The  CSDL  2  is  a  representative  LSS  for  a  precision  op¬ 
tical  space  system.  Eigenvalue  analysis  of  the  closed  loop  system  will  be 
used  to  evaluate  system  performance.  In  applying  the  control  method,  the 
output  (observation)  matrix  is  configured  to  supply  position  sensor,  rate 
sensor,  and  position/rate  sensor  information.  Point  force  actuators  pro¬ 
vide  the  state  variable  feedback  control. 

The  following  sections  will  detail  the  CSDL  2  model  and  its  finite  el¬ 
ement  representation.  The  control  and  matrix  transformation  methods  are 
discussed.  Finally,  the  computer  program  and  results  are  presented. 


II.  LSS  Model  Configuration 

A  central  issue  in  the  active  control  of  large  space  structures  is 
the  development  of  realistic  and  mathematically  correct  models  for  open 
and  closed  loop  dynamic  plants.  NASTRAN  and  SPAR  are  currently  the  pri¬ 
mary  tools  for  generating  models  of  sophisticated  conceptual  spacecraft. 

In  order  to  assess  the  control  technique  of  decentralized  direct  output 
feedback,  as  well  as  other  control  methods,  a  generic  model  is  required. 

The  Charles  Stark  Draper  Laboratory  (CSDL)  developed  just  such  a  model  as 
part  of  the  Active  Control  of  Space  Structures  (ACOSS)  program  sponsored 
by  DoD's  Advanced  Research  Projects  Agency  (4,5,14).  The  model  was  orig¬ 
inally  titled  ACOSS  Model  #2,  but  is  commonly  referred  to  as  the  CSDL  2 
model.  CSDL  2  is  a  finite  element  representation  of  an  optical  space 
structure.  It  has  incorporated  the  desired  features  of  structural  design 
based  on  realistic  sizes  and  weights,  a  simple  unclassified  optical  sys¬ 
tem  with  associated  performance  measures  and  tolerances,  and  a  set  of  dis¬ 
turbances  typical  of  equipment  vibration  and  attitude  control.  Figure  2 
shows  a  conceptual  view  of  the  structure  and  Figure  3  is  a  finite  element 
representation  of  the  model. 

CSDL  2  represents  a  three  mirror  optical  space  telescope  system. 

The  two  major  components  of  the  system  are  the  optical  support  structure 
and  the  equipment  section.  The  optical  support  structure  consists  of  the 
upper  mirror  support  truss,  the  lower  mirror  support  truss,  and  metering 
truss.  The  upper  mirror  support  truss  contains  the  primary  mirror  (convex 
surface)  and  the  tertiary  mirror  (concave  surface).  The  lower  mirror  sup¬ 
port  truss  contains  the  secondary  mirror  (flat  surface)  and  the  focal  plane 
(image  receiving  device).  The  metering  truss  maintains  the  mirror  separation 
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Figure  4.  CSDL  2  Optical  System 


and  is  the  key  section  when  examining  defocus.  The  optical  support  struc¬ 
ture  and  mirror  placement  are  shown  in  Figure  4.  Attached  to  the  lower 
side  of  the  lower  mirror  support  truss  is  the  equipment  section  which  con¬ 
sists  of  the  control  package,  modeled  as  a  rigid  body,  with  two  cantile¬ 
vered  flexible  solar  panels.  The  full  structure  is  approximately  twenty- 
eight  meters  high  and  has  a  mass  of  9300  kilograms.  The  structural  dimen¬ 
sions  are  shown  in  Figure  5. 

The  finite  element  model  of  the  structure  contains  fifty-nine  node 
points,  but  the  actual  structure  has  only  fifty-one  nodes.  The  extra 
nodes  were  added  to  provide  more  detail  in  the  modeling  of  the  mirrors  and 
equipment  section  (5).  The  coordinates  of  the  nodes  are  given  in  Table  I, 
and  the  placement  of  the  nodes  in  the  support  structure  are  shown  in  Fig¬ 
ure  6.  The  truss  members  are  fully  joined  so  that  bending  and  torsion  are 
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Figure  5.  Structural  Dimensions 


allowed.  The  truss  elements  are  made  of  graphite-epoxy  and  assumed  to  be 
massless.  The  system  mass  is  lumped  at  twenty-three  nodes  and  distributed 
as  shown  in  Table  II.  The  largest  mass  is  located  in  the  equipment  pack¬ 
age,  as  would  be  expected. 

This  model  makes  use  of  twenty-one  pairs  of  collocated  force  actua¬ 
tors  and  position  sensors.  A  list  of  sensor/actuator  locations  and  orien¬ 
tations  is  provided  in  Table  III. 

The  key  results  of  an  eigenvalue  analysis  performed  on  this  model 
are  listed  in  Table  IV.  The  generalized  mass  and  stiffness,  as  well  as 
the  natural  frequency  of  the  first  thirty-six  structural  modes  is  given. 

Of  the  modes  listed,  those  with  an  asterisk  are  considered  critical  for 

10 


-y-VL’i  •. 


,V_N  . '.  -•»  ■ 


.  .a  ~.  'ii  i 


Figure  6.  CSDL  2  Node  Placement 


LOS  performance  (15). 

Eigenvalue  analyses  will  provide  control  performance  information. 
Line-of-sight  performance,  as  well  as  defocus  along  the  z-axis,  are  impor 
tant  analysis  factors,  but  will  not  be  directly  addressed  in  the  results 
of  this  investigation.  Instead,  the  eigenvalue  analysis  will  provide  in¬ 
formation  on  controller  maintenance  of  modal  stability  and  on  controller 
Independence  (decoupling). 

The  controller  development  on  which  this  study  is  based  will  now  be 


presented. 
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TABLE  II 

CSDL  2  Lumped  Mass  Distribution 


67.4 

67.4 

67.4 

67.4 

69.5 
6.74 
69.5 
6.74 
6.74 
6.74 
69.5 
69.5 


44 

3500 

48 

81.9 

50 

1.638 

52 

73.8 

53 

73.8 

55 

163.8 

57 

81.9 

1001 

1000 

1002 

800 

1003 

1200 

1004 

600 

TABLE  III 

CSDL  2  Sensor/Actuator  Locations  and  Orientations 
(x,  y,  2  in  direction  cosines) 


Pair  Node  x  ^  z 

1  9  0  10 

2  9  0  0  1 

3  10  0  0  1 

4  11  1  o  0 

5  11  0  1  0 

6  11  0  0  1 

7  12  0  0  1 

8  27  1  0  0 

9  27  0  1  0 

10  27  0  0  1 

11  28  0  0  1 

12  29  0  1  0 

13  29  0  0  1 

14  30  0  0  1 

15  32  0  0  1 

16  33  0  0  1 

17  34  1  0  0 

18  34  0  1  0 

19  34  0  0  1 

20  35  0  1  0 

21  35  0  0  1 
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TABLE  IV 


Key  Results  of  NASTRAN 

Eigenvalue  Analysis  on  Nominal 

CSDL  2  Model 

Mode 

Generalized 

Mass 

Generalized 

Stiffness 

rad 
w  sec 

sec 

1-6* 

1.00E+00 

0.0 

0.0 

0.0 

7* 

1.00E+00 

5. 128E-01 

7. 161E-01 

5. 128E-01 

8 

1.00E+00 

8. 521E-01 

9.231E-01 

8.521E-01 

9 

1.00E+00 

8.835E-01 

9.399E-01 

8.835E-01 

10 

1.00E+00 

1.212E+00 

1. 101E+00 

1.212E+00 

11 

l.OOE+OO 

8.819E+00 

2.862E+00 

8. 189E+00 

12* 

l.OOE+OO 

1.266E+01 

3. 502E+00 

1.226E+01 

13* 

1.00E+Q0 

1.403E+01 

3.746E+00 

1.403E+01 

14 

1.00E+00 

1.492E+01 

3.863E+00 

1.492E+01 

15 

1.00E+00 

1. 599E+01 

3 . 998E+00 

1.599E+01 

16 

1.00E+00 

1.625E+01 

4 . 032E+00 

1.625E+01 

17* 

l.OOE+OO 

2.623E+01 

5. 122E+00 

2.623E+01 

18 

1.00E+00 

2.630E+01 

5. 128E+00 

2.630E+01 

19 

1.00E+00 

2.677E+01 

5. 174E+00 

2.677E+01 

20 

1.00E+00 

3.310E+01 

5.753E+00 

3.310E+01 

21* 

1.00E+00 

3.730E+01 

6. 197E+00 

3.730E+01 

22* 

l.OOE+OO 

5.301E+01 

7.281E+00 

5.301E+01 

23 

1.00E+00 

9.498E+01 

9.746E+00 

9.498E+01 

24* 

l.OOE+OO 

1.241E+02 

1. 114E+01 

1.241E+02 

25 

l.OOE+OO 

1. 999E+02 

1.414E+01 

1 . 999E+02 

26 

l.OOE+OO 

2.001E+02 

1.416E+01 

2.001E+02 

27 

l.OOE+OO 

4. 654E+02 

2. 157E+01 

4.654E+02 

28* 

l.OOE+OO 

4.705E+02 

2. 169E+01 

4.705E+02 

29 

l.OOE+OO 

6. 182E+02 

2.468E+01 

6. 182E+02 

30* 

l.OOE+OO 

6.275E+02 

2. 505E+01 

6.275E+02 

31 

l.OOE+OO 

6.481E+02 

2 . 546E+01 

6.481E+02 

32 

l.OOE+OO 

7.428E+02 

2.725E+01 

7 . 428E+02 

33 

l.OOE+OO 

1.700E+03 

4. 123E+01 

1.700E+03 

34 

l.OOE+OO 

2.568E+03 

5.067E+01 

2 . 568E+03 

35 

l.OOE+OO 

2.821E+03 

5.311E+01 

2.821E+03 

36 

l.OOE+OO 

3.095E+03 

5. 563E+01 

3.095E+03 
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III.  System  Model 


Equations  of  Motion 

As  presented  by  Calico  and  Miller  (2,13),  the  system  model  may  be 
developed  from  the  vibrational  equations  of  motion  for  a  large  space  struc¬ 
ture  given  generally  as 

Mg  +  Eg  +  Kg  =  Du  (1) 

where  M  is  a  nxn  symmetric  mass  matrix,  E  is  a  nxn  sytwnetric  damping  ma¬ 
trix,  K  is  a  nxn  symmetric  stiffness  matrix,  D  is  a  nxm  matrix  of  nodal, 
attitude-evaluated  actuator  locations,  g  is  a  nxl  generalized  coordinate 
vector,  and  u  is  a  mxl  control  input  vector.  Introducing  the  nxn  modal 
matrix  $  for  Eq  (1),  such  that 

9  *  $n  (2) 

where  ri  is  the  n-vector  of  modal  coordinates,  Eq  (1)  may  be  written  as 

^ij  n  +  n  +  n  =  ^Du  (3) 

the  a),  and  terms  being  natural  frequencies  and  damping  coefficients, 
respectively,  of  the  specific  modes.  The  properties  of  the  modal  matrix 
$  are  such  that  the  coefficients  of  Eq  (3)  are  given  be 

[  I  ]  =  $tm  $ 

]  =  $TE  $  (4) 

^  oj2  J  s  $ 
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_I  =  nxn  identity  matrix 

2eui  =  nxn  diagonal  damping  matrix 

co2  =  nxn  diagonal  matrix  of  eigenvalues  of  Eq  (1) 

Equation  (3)  may  now  be  converted  into  a  state  space  representation 
of  the  system,  given  by 

x  =  Ax  +  Bu  (5) 

in  which  A  is  the  nxn  plant  matrix,  B  is  the  nxm  input  matrix,  x  is  the 
nxl  state  vector,  and  u  is  the  mxl  control  input  vector.  These  system 
parameters  in  state  space  are 


The  complete  state,  however,  is  normally  not  available,  so  Eq  (5)  must  be 
supplemented  by  an  output  equation.  State  space  form  gives  the  sensor 
output  as 


when  both  position  (p  subscript)  and  velocity  (v  subscript)  sensors  are 
used.  Expressing  Eq  (7)  with  the  state  vector  x 

y  =  c*  (8) 

where 


Equations  (5)  and  (8)  form  the  large  space  structure  model  available  to 
the  control  designer.  These  equations  will  be  further  explained  so  they 
will  hold  more  significance  when  being  applied  to  modal  control  of  flexi¬ 
ble  structures. 

Control  Model 

The  full  structural  model  is  represented  by  the  2n-dimensional  state 
vector  x.  As  noted  earlier,  it  is  impossible  to  model  all  of  the  possible 
modes  for  a  complex  structure,  and  of  those  modeled,  even  fewer  will  be 
actively  controlled.  Assuming  that  multiple  controllers  are  available, 
each  controlling  a  small  subset  n.  of  nodes,  as  in  this  investigation,  the 
state  vector  may  be  simply  represented  by 

x  =  jxl,  xl . .  xj,  xj,  xjm  J  T  (10) 

The  x^  terms  represent  2n. -vectors  of  modal  amplitudes  and  velocities  as 
defined  by  the  last  of  Eq  (6)  controlled  by  the  i  th  controller  of  N  con¬ 
trollers  present.  The  xr  represents  a  2nr-vector  of  modeled  residual  modes 

and  the  x..„  represents  a  2n  -vector  of  unmodeled  modes, 
urn  r  urn 

As  defined  earlier,  the  unmodeled  modes  are  those  which  exist  but 
are  beyond  the  number  of  modes  in  the  structural  model.  These  will  no 


longer  appear  in  the  derivations.  The  residual  modes  are  those  which  are 
modeled  but  not  controlled.  The  controlled  modes  are  those  which  require 
active  control  in  order  to  obtain  satisfactory  system  response.  The  se¬ 
lection  of  the  modes  to  be  controlled  and  their  assignment  to  one  of  the 
N  controllers  is  left  to  the  control  designer. 

It  should  be  noted  that  so  far  suppressed  modes  have  not  directly 
been  referred  to  even  though  they  were  defined  earlier  in  the  text.  They 
have  not  been  ignored,  in  fact  they  are  included  within  the  controlled  modes, 
in  the  following  manner:  In  a  multiple  controller  design,  the  individual 
controller  actively  controls  those  modes  assigned  to  it.  But  the  modes 
assigned  to  the  other  controllers  interact  with  this  individual  controller, 
causing  control  and/or  observation  spillover.  In  the  process  of  controlling 
the  system,  each  controller  may  contribute  to  the  elimination  of  observation 
and  control  spillover  in  the  system.  Thus  each  controller,  in  effect  "sup¬ 
presses"  the  modes  contained  in  the  other  controllers.  In  other  words, 
the  controlled  modes  of  one  controller  are  the  suppressed  modes  of  another 
controller.  Therefore,  the  suppressed  modes  are  contained  implicitly  within 
the  controlled  modes.  So,  like  the  unmodeled  modes,  the  suppressed  modes 
exist  in  the  system,  but  will  not  be  mentioned  in  the  derivations  since  they 
are  included  implicitly  in  the  controlled  modes. 

Continuing  with  the  derivation,  the  notation  of  Eq  (10)  may  be  used 
to  express  the  state  equations  as  follows 

xi  a  A..X.J  +  B.u,  i  =  1,  2 . .  N  (11) 

xr  =  Arxr  +  Bfu  (12) 
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N 

y  =  2 

i= 


where  the  A,  B,  and  C  matrices  are 


j  =  1,  2,  . . . ,  N,  r  (14) 


j  =  1(  2 . N,  r  (15) 


cj  =  [Cpj$  :  Cvj$  ]  j  =  1.  2,  ....  N,  r  (16) 

Moreover,  the  lower  partition  of  Eq  (15),  the  4TD.j  and  4TDp  matrices,  are 


of  the  form 


where 


4  D.  = 

1  l 


$u„  =  4L 
r  r 


-  h  ■  aik 
<Vjk  -  h  ■  3rk 


The  4j  are  the  column  vectors  of  the  matrix  4  and  the  d.^  and  drk  are  the 
column  vectors  of  the  and  Dp  matrices,  respectively.  Using  the  forms 
given  in  Eqs  (17)  to  (20),  the  4TD.  and  4TDp  matrices  may  be  represented 


Vi 


)  1 1 

(^.j)l2 

a 

V.  = 

1 

(^^ ) 2  X 

(^^22 

a 

(21) 

(rt.) 

wvn.i 

(♦,). 

^  n^2 

^i^n.n 

i  a 

(^r)u 

(^r)l2 

•••  (*r>in 

a 

¥  = 

r 

(^r)21 

(^r)22 

($  ) 
VHr'2n, 

a 

(22) 

(tf,  ) 

(ij,  ) 

'vr#nr2 

-  (*>V. 

J 

where  n.  is  the  number  of  inodes  in  the  i  th  controller,  n  is  the  number 
1  r 

of  residual  modes,  and  nfl  is  the  number  of  actuators  employed.  This  al¬ 
lows  Eq  (15)  to  be  rewritten  as 


j  =  1,  2 . N,  r 


(23) 


In  simpler  terms,  the  rows  of  the  matrices  represent  the  amplitude  of 

w 

each  structural  mode  along  the  line  of  action  of  each  actuator  location. 

The  dimension  of  the  matrices  is  n.  x  n.  making  the  di- 

J  J  a 

mension  of  the  B.  matrices  2n.  x  n.  when  the  upper  null  partition  is  in- 
j  j  a 

eluded.  Likewise  it  can  be  seen  that  the  C  and  C  .$  partitions  of  Eq 

PJ  vj 

(16)  are  of  dimension  n$  x  n^.  where  ng  is  the  number  of  sensors  employed. 
This  makes  the  dimension  of  the  C.  matrices  n  x  2n.. 

J  ^  J 


Examining  the  C.  matrix  more  closely,  the  C  .  and  C  .  terms  are  the 

J  r J  *  J 


position  and  velocity  coefficient  matrices,  respectively,  of  the  sensors 
employed,  assuming  that  both  position  and  velocity  sensors  are  used.  How¬ 


ever,  if  only  rate  sensors  are  used,  this  makes  the  C  .  into  zero  matrices 

pj 


so  that  Eq  (16)  now  becomes 

Cj  =  [°  :  Cvj*]  j  =  u  2 . N’ 


(24) 


Furthermore,  when  collocated  sensor/actuators  are  employed  with  the  same 
alignment,  this  simplifies  even  more.  In  this  special  case 


so  that 


c  .*  = 

vj* 

[•’*,  ]  ’ 

(25) 

C  •$  = 
vr 

>j] 

(26) 

Vi 

j  *  1,  2,  ...,  N,  r 

(27) 

This  simplicity  is  the  prime  advantage  of  using  either  rate  or  position 
sensors  only.  As  pointed  out  in  the  B.  matrix,  the  columns  of  the  v  •  ma- 

J  J 

trix  in  C.  represent  the  rates  of  each  structural  mode  at  each  sensor  lo- 
cation  along  the  line  of  the  sensor. 

The  equations  thus  derived  are  very  general  in  form  and  are  indepen¬ 
dent  of  structural  complexity.  Only  the  matrix  dimensions  will  vary  de¬ 
pending  on  the  number  of  sensors,  actuators,  and  modes  studied.  This  gen¬ 
eral  development  can  lead  one  to  understand  the  wide  variety  of  structures 
to  which  the  following  analysis  may  be  applied. 


Modal  Control 


The  controller  design  for  N  controllers  will  be  based  upon  the  model 
given  by  Eqs  (11),  (12),  and  (13)  which,  restated,  are 


=  A..X..  +  B.u 

(ID 

=  +  B  u 

r  r  r 

(12) 

N 

=  ACiXi  +  Crxr 

(13) 

The  term  Bru  in  Eq  (12)  represents  control  spillover  from  the  residual  con¬ 
trol  matrix.  Similarly,  the  term  Crxf  is  called  observation  spillover. 

These  spillovers  can  destabilize  the  residual  modes.  In  decentralized  con¬ 
trol  systems  another  set  of  spillovers  appears.  If  the  states  are  available 
to  the  control  designer,  the  feedback  controller  may  be  designed  such  that 

N  „ 

u  =  Z  K.x.  (28) 

i=i  1  1 

A  A 

where  K.  are  the  control  gain  matrices.  The  development  of  the  K.  matrices 
will  be  presented  shortly.  Now  the  residual  state  equation  becomes 


N  „ 

B  Ar*r  +  1*1 


(29) 


which  is  spillover  from  the  controlled  modes  to  the  residuals.  The  con¬ 
trolled  s.tate  equation  is 


Xi  =  A..X.  +  l  KjXj 


(30) 


and  shows  the  control  spillover  from  controller  to  controller. 


22 


Now  for  direct  output  feedback  the  control  vector  is 


(31) 

(32) 

(33) 


and  in  matrix  form  for  the  N  controller  with  residuals 

Ai+BjkCi  B!<C2  •••  Bt<CN  B^C^ 

B2kCi  A2+B2<C2  •••  B2icC|y  B2<Cr 

•  «  •  •  x  (34) 

•  •  •  • 

Bn<Cx  BnkC2  •••  ^+BnkCn  BN<Cp 

Br<Ci  Br<C2  •••  Br<CN  Ar+Br<cr 

N  „ 

where  k  =  E  (K. ).  The  last  column  in  this  plant  matrix  represents  obser- 
i=i  1 

vation  spillover;  the  last  row,  control  spillover.  The  upper  off-diagonal 
terms  are  the  controller  to  controller  observation  spillover.  The  lower 
off-diagonal  terms  the  controller  to  controller  control  spillover  terms. 

To  assure  stability  of  the  system  based  on  the  stability  in  each  controller, 
the  system  matrix  must  be  within  a  similarity  transformation  of  either  upper 
or  lower  block  triangular.  Lower  block  triangular  form  requires 

B.K,  =  0  and  K.C.  =0  i  =  1,2,...,N-1;  j  =  i  +  1 . N  (35) 


and  upper  block  triangular  form  requires 


B.Kj  =  0  and  K.C.  -  0  i  =  j  +  j  =  1,2,...,N-1  (36) 

These  conditions  are  identical  to  those  required  by  the  decentralized  op¬ 
timal  controller  as  found  by  Calico  and  Miller  (2). 

The  block  structure  of  Eq  (34)  can  be  more  easily  seen  in  a  specific 
example.  The  three  controller  case  will  now  be  examined. 

Three  Controllers 

Setting  N  equal  to  3  and  following  the  form  given  in  the  previous 
development  for  N  controllers,  the  state  equations  for  a  three  controller 
system  are  given  as 


Xi  =  AxXx  +  BiU 

(37) 

• 

x2  *  A2x2  +  B2u 

(38) 

• 

x3  =  A3X3  +  B3u 

(39) 

• 

Xr  *  Arxr  +  Bru 

(40) 

The  control  applied  is 


u  =  z  K,y  =  ZKAZ  C  x  +  C  x  )  (41) 

i=i  1  i=i  1  m=i  m  m  r  r 

u  =  (Kj  +  K2  +  KsMCxXx  +  C2x2  +  C3x3  +  Crxr)  (42) 

Substituting  Eq  (42)  into  (37)  through  (40),  the  system  equations  are 


x,  *  AiXx  +  Bj(Ki  +  K2  +  K3 ) ( C2 x !+  C2x2+  C3x3+  Crxr) 


1 


4 


s  (Aj  +  Bj<Ci )xi  +  Bx<C2x2  +  Bj^C3X3  +  Bi<C  x, 


(43) 


(44) 


x  2  =  A2X2  +  B2(Ki  +  K2  +  K3HC1X1  +  C2x2  +  C3X3  +  Crxr ) 

=  B2kC1Xi  +  ( A2  +  B2kC2)x2+  B2^C3x3  +  B2<Cr>xr 

x3  =  A3X3  +  B  ( Kx  +  K2  +  K3)(CiX!  +  C2x2  +  C3X3  +  Crxr) 

=  B3icCiXi  +  B3tcC2X2  +  ( A3  +  B3kC3)x3  +  B3KCrxr  (45) 

xr  *  Afxf  +  Bf(Ki  +  K2  +  K3)(C1xl  +  C2x2  +  C3X3  +  Crxr) 

*  B^CxX!  +  BrKC2X2  +  BrKC3x3  +  (Ar  +  Br<Cr)xr  (46) 


Now  combining  these  system  equations  with  an  augmented  state  vector  x 


=  {  XjT,  x2T,  x3T,  xrT  |  T 


the  closed  loop  system  equation  may  be  giver  a* 


(47) 


x  = 


Ai  +  Bi<Ci 

B2kCi 

B3kCi 


Br<Cl 


B1KC2 

a2  +  b2<c2 

B3KC2 


Bi<C3 

B2kC3 

A3  +  B3KC3 


Br<C3 


Bi<C, 


B2^C„ 


B3KC 


A  +  B  <C 
r  r  r 


(48) 


Inspection  of  Eq  (48)  shows  that  individual  controller  stability  cannot 
guarantee  overall  system  stability.  Block  triangularization  may  be  used 
to  remove  controller  to  controller  spillover,  but  the  spillover  from  the 
residuals  will  remain.  Triangularization  requires  the  satisfaction  of  Eq 
(35)  conditions.  These  conditions  may  partially  be  met  by  using  selective 
sensor  and  actuator  placements,  but  in  general  their  satisfaction  requires 
careful  consideration. 
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The  approach,  used  In  this  study,  proposes  that  transformation  matrices 


T.  and  r.  be  found  such  that  B.T.  =  0  and  r.C.  =  0  ,  j>i  for  lower 

J  '  *  J  *  J 

block  triangular  form,  j<i  for  upper  block  triangular  form.  It  will  be 
shown  shortly  how  the  T.  and  r.  may  be  used  in  eliminating  destabilizing 

J  * 

spillovers. 

Design  Feedback  Gain 

The  suboptimal  output  feedback  gains  for  this  study  are  obtained  using 
an  extension  of  the  Kosut  approximation  (18).  The  state,  output,  and  con¬ 
trol  equations  are 

x  =  Ax  +  Bu  (49) 

y  =  Cx  (50) 

u  =  Ky  =  KCx  (51) 

The  direct  output  feedback  allows  us  to  rewrite  Eq  (49)  as 

x  =  (A  +  BKC)x  ’  (52) 

A 

It  is  desirable  for  the  closed  loop  system  matrix,  A  +  BKC,  to  equal  the 
closed  loop  system  matrix  whose  control  gain  matrix  is  the  optimal  steady 
full-state  gain,  denoted  as  G.  This  gain,  G,  is  determined  using  optimal 
control  (3).  The  control  gain  matrix  G  is  derived  by  defining  a  quadratic 
performance  index  J  such  that 


oo 

J  ~  h  J  (xTQx  +  uTRu)dt 


where  Q  is  an  nxn  positive  semidefinite  weighting  matrix  and  R  is  an  mxm 
positive  semidefinite  weighting  matrix.  It  is  desired  to  minimize  this  in¬ 
dex  subject  to  Eq  (11).  Then  the  optimal  solution  to  the  minimization 


problem  is 


G  =  -R_1BTS  (54) 

where  S  is  the  solution  to  the  steady  state  matrix  Riccati  equation 

SA  +  ATS  -  SBR'XBTS  +  Q  =  0  (55) 

The  control  applied  with  optimal  feedback  gains  becomes 

u  =  Gx  (56) 

for  which  the  closed  loop  system  is 

x  =  (A  +  BG)x  (57) 

/\ 

The  two  systems  are  equated  and  the  suboptimal  gain  matrix  K  is  solved  for 

(A  +  BKC)x  *  (A  +  BG)x  (58) 

A  +  BKC  =  A  +  BG  (59) 

BKC  =  BG  (60) 

KC  =  G  (61) 


In  general,  the  output  matrix  is  not  square  and  at  this  point  the  Kosut 

A 

approximation  solves  for  K  using  a  pseudo-inverse  method 


~  T  T 

KCC  =  GC1 

(62) 

k(cct)(cct)‘1  =  gct(cct)-1 

(63) 

K  =  GCT(CCT)'1 

(64) 

This  solution  is  valid  if  the  output  matrix  is  of  full  rank  making  CCT  non- 

/V 

singular  and  invertible.  For  the  case  where  C  is  not  of  full  rank,  K  can 


•V 
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still  be  found  by  using  the  inverse  of  the  singular  value  decomposition  of 
C,  denoted  as  C+  (22).  This  inverse  is  also  known  as  the  generalized  or 
Penrose  inverse.  The  method  of  singular  value  decomposition  is  discussed 
in  the  next  chapter.  Armed  with  this  generalized  inverse  and  solving  for 
K  from  Eq  (61) 


A  X  X 

KCC  =  GC 

~  X 

K  =  GC 


(65) 

(66) 


This  relatively  simple  procedure  allows  us  to  determine  the  suboptimal  gains 
for  an  output  matrix  of  any  rank. 

These  equations  can  be  generalized  for  decentralized  controllers  and 
allow  us  to  complete  the  development  of  the  transformed  DOFB. 


Transformed  DOFB 

Consider  an  output  transformed  by  I\ 


and  a  control  law 


w.  =  r.y. 
i  r  i 


N 

u  =  Z  T.v. 
i=i  1  1 


(67) 


(68) 


where  =  K*w..  When  r.  and  T.  are  such  that  the  conditions  B.jT.  =  0 
and  r.C.  =  0  for  j>i  or  i<j  are  met,  destabilizing  controller  spillover 

*  J 

is  eliminated.  Either  of  these  conditions  will  lead  to  a  closed  loop  state 
system  where  the  controller  matrices  on  the  diagonal  are  decoupled  and  are 
the  result  of  solving  the  following  open  control  problems  (6): 


x.  =  Ax.  +  B?v. 


(69) 


(70) 


w. 

1 


B? 


BiTi 


c* 


r  c 

iLi 


vi 


Ki*i 


(71) 

(72) 

(73) 


The  gain  K*  is  obtained  by  solving  for  a  transformed  optimal  gain  G|.  G| 
is  found  by  substituting  into  Eqs  (54)  and  (55)  and  finding  the  gener¬ 
alized  inverse  of  C*.  Hence,  the  transformed  output  gain  is 


K*  =  G*C* 

the  design  feedback  control  gain  matrix  K  is  given  by 


K,-  *  Wi 


and  the  closed  loop  equation  becomes 


N  _ 

=  A..  +  B.(  E  K..  )C. 


In  matrix  form  Eq  (76)  is 


A i +B i K  x  C i 

0 

0 

•  •  •  0 

BjACr 

B2K1C1 

A2+B2K2C2 

0 

•  •  •  0 

B2ACr 

B3(Kl+K2)Ci 

B3K2C2 

• 

• 

• 

fl 

• 

• 

• 

• 

• 

• 

• 

* 

0 

• 

N-i„ 

Bn(  2  K. )Ci 

"  i  =  i  1 

N-2, 

Bn(  2  K.)C2 

n  i=i  1 

•  •  • 

am+bwLcm 

N  N  N  N 

VCr 

BfACi 

BrAC2 

•  •  ■ 

B  ACm 
r  N 

Ar+BrAC. 

29 


-1-—  ■  - 1-  •  -Ti  - - - ^.L- 


(74) 


(75) 


(76) 


(77) 


where  A  =  l  K. .  For  a  three  controller  system  with  residuals,  which  this 

i  =  i  1 

study  implements,  Eq  (77)  becomes 


Ai+BiKiCi 

0 

0 

BiACr 

B2K1C1 

A2‘*"82  K2C2 

0 

B2ACp 

X 

B3(Ki+K2)C1 

B3K2C2 

A3+B3  K3  C3 

B3ACr 

X 

B/Ci 

BrAC2 

BrAC3 

WCr_ 

A  four  controller  system  without  residuals  is  also 

examined 

and 

loop  system  equation  with  design  gains 

should  look 

1  ike 

Ai  +Bi  Ki  Ci 

0 

0 

0 

BzKiCi 

A2+B2 K2  C2 

0 

0 

x  = 

B3 (K1+K2 )Ci 

83  K2  C2 

A3+B3  K3  C3 

0 

B4(K1+K2+K3)C1 

B4(K2+K3)C 

2  B4K3C3 

A4+B4K4C4  1 

(78) 


*  (79) 


It  is  quite  obvious  from  comparing  Eqs  (78)  and  (79)  that  the  residual  terms 
can  have  a  significant  effect  on  system  eigenvalues  and  hence  system  stabil¬ 
ity,  while  the  system  eigenvalues  of  Eq  (79)  will  always  be  those  of  the  in¬ 
dividual  controllers  along  the  diagonal.  With  the  control  developed,  we  need 
to  establish  how  many  modes  may  be  controlled  in  a  given  controller,  i.e. 
define  the  sensor/actuator  requirements. 


Sensor/Actuator  Requirements 

As  mentioned  for  three  controllers,  in  order  to  perform  spillover  sup 
presslon,  one  or  more  gain  matrices  must  be  made  orthogonal  to  N-l  B  or  C 

A 

matrices.  For  example,  from  Eq  (35),  to  satisfy  the  expression  for  B.K. 


the  columns  of  Ki  must  be  simultaneously  orthogonal  to  the  rows  of  B2  through 
B^.  In  other  words,  the  columns  of  Kx  must  be  in  the  null  space  of  the  ma¬ 
trix  B2N,  where  B2N  is  defined  as 


The  null  space  of  B«w  has  dimension  P9W  given  as 


P2N  ^na  "  r2fP 


where 


na  =  number  of  actuators 

cl 


r2N  =  rank  of  B2N  <  min  (n2  +  n3  +  ...  +  n^,  nfl) 


Therefore,  Ki  has  P2N  columns. 

The  number  of  actuators  must  exceed  the  rank  of  B2N  in  order  for 

A 

B2NKl  =  <>the™ise  the  system  is  overspecified  and  no  transformation 

A 

matrix  exists  which  will  drive  the  B^Ki  to  zero.  If  rows  of  B2N  are  lin¬ 
early  independent,  then  the  number  of  actuators  needed  is  given  as 


na  >  E  n. 
a  1=2  1 


and  if  the  rows  are  not  linearly  independent,  na>r2N<  It  can  be  seen  that 
the  other  control  gain  matrices  will  have  a  sufficient  number  of  actuators 
if  the  Inequality  in  Eq  (82)  is  met.  A  similar  study  shows  that  the  number 


of  sensors  needed  is,  for  CgN  with  linearly  independent  columns 

N-i 

n.  >  in.  (83) 

s  i=i  1 

and  for  columns  that  are  not  linearly  independent,  ng  >  r^  here  is  the 
rank  of  Cgj^- 

Likewise,  for  the  lower  block  triangular  conditions  given  in  Eq  (51), 
and  actuator  and  sensor  requirements  can  be  shown  to  be 


N-i 

(84) 

N 

Z  n 

1  =  2 

(85) 

for  linearly  independent  rows  and  columns  of  the  B2N  and  C2N  matrices,  re¬ 
spectively. 

It  should  be  pointed  out  that  satisfying  the  inequalities  given  in 
Eqs  (82)  through  (85)  may  actually  require  more  actuators  and  sensors  than 
indicated.  As  an  example,  consider  a  thirty  mode  model.  Using  three  con¬ 
trollers,  each  with  ten  modes,  the  above  inequalities  require  at  least  twenty 
sensors  and  twenty  actuators  to  assure  decoupled  system  stability.  However, 
to  control  and  observe  the  system,  at  least  one  more  sensor  and  actuator 
would  be  required.  These  conditions  must  be  met  in  order  to  generate  the 
transformation  matrices  and  to  implement  the  controllers.  With  these  re¬ 
quirements  and  control  model  developed,  the  technique  to  obtain  the  trans¬ 
formation  matrices  and  are  presented. 


IV.  Transformation  Technique 


The  transformation  technique  used  in  this  study  is  the  same  one  ap¬ 
plied  by  Aldridge  and  Miller  (1,13).  It  has  been  mentioned  several  times 
that  the  closed  loop  state  equations,  Eqs  (11)  and  (33),  will  be  put  into 
block  triangular  form  by  the  selective  elimination  of  control  and  obser¬ 
vation  spillover  terms.  However,  the  exact  details  of  this  spillover  elim 
ination  have  been  neglected  until  now.  The  following  will  describe  the 
generation  of  the  transformation  matrices,  which  were  referred  to  speci¬ 
fically  as  T  and  r  in  the  previous  section.  The  T  matrix  is  a  transfor¬ 
mation  matrix  for  the  control  spillover  and  the  r  matrix  is  a  transfor¬ 
mation  matrix  for  the  observation  spillover. 

In  a  single  controller  case,  it  can  be  seen  that  the  spillover  terms 
which,  when  eliminated,  will  assure  system  stability  are  B$K  or  KCS,  where 
the  s  subscript  designates  modes  to  be  suppressed.  An  immediately  obvi¬ 
ous  solution  to  this  is  K  =  0.  But  this  solution  will  make  the  respec¬ 
tive  terms  B  K  and  KC  equal  to  zero  also.  Therefore,  this  solution  is 
unsatisfactory.  The  transformation  method  generates  a  solution  which  for 
a  single  controller,  is  subject  to  the  conditions 


BSK  =  0 
KC$  =  0 


while  maintaining 


[»j 


It  would  also  be  desirable  to  apply  the  additional  constraint  to  the  re¬ 
sidual  modes: 


BrK  =  0 

(90) 

KCf  =  0 

(91) 

however,  due  to  the  large  number  of  structural  modes  present  in  the  model, 
this  constraint  is  not  realistic  and  will  be  ignored  in  this  development. 
The  effects  of  the  residual  spillovers  may  be  minimized  by  the  careful  se¬ 
lection  of  modes  designated  as  residual  or  suppressed,  so  as  to  create  a 
frequency  separation  between  the  residuals  and  the  bandwidth  of  the  con¬ 
troller. 

For  a  multiple  controller,  the  conditions  given  in  Eqs  (86)  and  (87) 
apply,  but  now  the  B  and  C  matrices  may  take  on  the  form  illustrated  by 
Eq  (80)  and  will  be  referred  to  as  the  B.^  and  C-N  matrices.  Instead  of 
discussing  all  of  the  possible  combinations  of  B^N  and  C^N  for  N  control¬ 
lers,  take  as  an  example  the  first  condition  given  in  Eq  (35),  that  is 

B.K.  =0  j  =  1,2,... ,N-1;  i  =  j+l,...,N  (92) 

Given  N  controllers,  Kf  will  have  to  be  made  orthogonal  to  N-l  B^  matrices. 
The  N-l  B.  matrices  may  be  combined  into  a  single  matrix  such  that 


Therefore,  one  of  the  conditions  to  be  met  is  B^Kx  =  0.  In  other  words, 
the  Ki  matrix  must  be  transformed  such  that  its  columns  are  orthogonal  to 
the  rows  of  B2N,  or  as  seen  in  the  previous  section,  Kx  must  be  in  the 
null  space  of  B2^.  This  is  the  most  difficult  case  for  N  controllers  and 
is  presented  only  to  describe  how  the  multiple  matrix  is  set  up.  The  re¬ 
mainder  of  the  derivation  will  be  in  terms  of  a  generic  B$  matrix  which 
represents  B2N  and  B2  alike. 

The  transformation  matrix  sought  will  be  referred  to  as  T  and  will 
be  such  that 


BST  =  0 


B$  has  the  row  dimension  of  nm  (the  number  of  modes  to  be  suppressed)  and 
the  column  dimension  of  n.  (the  number  of  actuators).  T,  therefore,  has 
dimensions  of  n  by  na  -  n.  If  there  are  fewer  actuators  than  linearly 

a  a  in 

independent  modes,  then  no  solution  matrix  T  exists.  The  system  is  over- 
specified  in  this  case,  meaning  there  are  more  equations  than  unknowns. 

If  the  number  of  modes  and  actuators  are  equivalent,  then  the  system  is 
stable  but  uncontrollable  (or  in  the  case  of  the  KC  matrices,  the  system 
is  stable  but  unobservable  assuming  an  equal  number  of  sensors  and  actu¬ 
ators).  The  actuators  (sensors)  are  saturated  with  maintaining  stability 
alone.  Simply  stated,  the  conditions  given  in  Eq  (82)  and  (83)  must  be 
met  in  order  to  generate  a  transformation  matrix. 

To  illustrate  the  result  of  applying  the  transformation  matrix  de¬ 
scribed  above,  consider  the  following  system  of  controlled  and  suppressed 
modes: 


x.  =  A  x„  +  B  u 
c  c  c  c 


35 


(96) 


where 


x.  * 


Ax  + 

Vs 


V 


u  = 


Mr 
c  c 


(97) 


The  Bsu  term  of  Eq  (91)  is  a  control  spillover  term  which  may  be 
adversely  affected  by  the  control  applied  to  Eq  (96).  The  elimination  of 
this  term  requires  the  use  of  a  transformation  matrix  T  such  that 

BST  =  0  (98) 

while  maintaining 

BCT  f  0  (99) 


This  transformation  matrix  may  be  used  to  define  a  new  control  v  as 

u  =  Tv  (100) 

Inserting  this  expression  into  the  state  equations  given  in  Eqs  (91) 
and  (92)  yields 

*c  -  Vc  +  V*  (101) 

K  -  Vs +  v5  <102> 


Letting  BCT  =  B*  and  knowing  that  BgT  =  0  the  new  system  is  described 
by 


-  Vc +  V 

(103) 

=  Vs 

(104) 

in  which  no  controller  spillover  exists.  The  new  control  vector  will  be 
shown  to  be 


(105) 


With  this  general  overview  of  the  transformation  process  and  its  re¬ 
sults,  the  development  of  the  matrix  T  will  now  be  discussed. 

The  major  tool  used  obtaining  this  result  is  called  the  Singular 
Value  Decomposition  (22).  The  matrix  to  be  decomposed  is  B$  which  has  di¬ 
mensions  n,  x  n,  and  can  be  described  by 
m  a 

Bs  =  WEVT  (106) 

where  W  is  an  (nm  x  nm)  orthogonal  matrix  of  left  singular  vectors,  V  is 
an  (n,  x  n_)  orthogonal  matrix  of  right  singular  vectors,  and  E  is  an 

a  a 

(nffl  x  na)  matrix  with  the  s  singular  values  of  B$  in  the  first  s  entries 
along  the  main  diagonal  and  zeroes  in  all  other  positions: 


The  total  number  of  singular  values  present  is  equal  to  the  rank  of 

the  Bg  matrix,  and  they  are  all  non-negative.  Assuming  B$  is  of  full  rank, 

the  dimension  of  S  *  min  (n  ,  n). 

a  m 


By  partitioning,  the  W  matrix  can  be  defined  by 


M  = 


(109) 


where  W$  is  an  n^  x  s  matrix  of  left  singular  vectors  associated  with  the 
non-zero  singular  values,  Wr  is  an  nffl  x  r  matrix  of  left  singular  vectors 
associated  with  the  zero  singular  values,  and 


s  +  r  =  nm 
m 


(110) 


Similarly,  by  partitioning,  the  V  matrix  can  be  defined  by 

v  *  [».;  VP]  (U1) 
where  V_  is  an  s  x  na  matrix  of  right  singular  vectors  associated  with  the 

5  a 

non-zero  singular  values,  V„  is  a  p  x  n.  matrix  of  right  singular  vectors 

p  a 

associated  with  the  zero  singular  values,  and 

s  +  p  =  na  (112) 


Remembering  the  V  matrix  is  orthogonal  and  noting  then  that 

VsTVp  =  0  (113) 

the  decomposed  matrix  may  be  written  as 

bsT  ■  wss  vsT  vp  ■  BS  \  -  0  <u4> 

which  leads  to  the  conclusion  that  the  transformation  matrix  desired  is 
given  In  the  matrix  of  right  singular  values  associated  with  the  zero 
singular  values: 


(115) 


T  *  V. 


'  N 


where  T  f  0. 

Once  the  transformation  matrix  is  found,  implementation  is  relatively 
simple.  Equation  (100)  defined  a  as  the  new  control  input  vector. 

This  is  now  given  as 


v,  =  G*  X, 


(116) 


The  Gt  matrix  is  found  in  the  same  manner  as  in  the  Modal  Control  subsec¬ 
tion  of  Section  III  in  which 


where 


Gi 


RT 


T1Tr1T1 


R..  is  the  positive  definite  weighting  matrix  in  Eq  (53) 


(117) 


(118) 


-  B,!, 


(119) 


and  is  the  solution  to  the  matrix  Riccati  equation 


siN ♦  AiT$i  -  siB?Rr'B?Tsi +  ■  ° 


(120) 


Simple  manipulation  of  Eq  (117)  will  show  that  the  transformed  gain  ma¬ 
trix  is  finally  given  by 


G*  ■  T.G. 

i  ii 


(121) 


■V- 

*  V„ 


Substitution  of  this  back  into  the  state  equations  yields  a  closed  loop 
system  which  is  block  triangular  with  no  control  spillover. 
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This  technique  may  be  paralleled  to  obtain  a  r  transformation  matrix 
for  the  observation  gain  to  eliminate  observation  spillover.  Substituting 
C$T  for  Bs>  and  r  for  T  will  give  the  same  results  with  r  equal  to  V  . 
Otherwise,  a  new  derivation  using  C$  will  give  the  result  that  the  obser¬ 
vation  transformation  matrix  r  is  equal  to  the  transpose  of  the  matrix  Wr 
of  left  singular  vectors  associated  with  the  zero  singular  values  of  Cs. 
Here,  as  alluded  to  earlier,  the  number  of  sensors  must  be  greater  than 
the  number  of  modes  to  be  suppressed. 


V.  Computer  Model 


Simplicity  and  flexibility  were  considered  in  the  development  of  the 
computer  program.  The  general  format  given  by  Aldridge  (1)  was  followed 
as  it  provided  a  simple  progression  of  logic.  The  program  was  modified 
to  calculate  the  suboptimal  gains  and  accommodate  a  design  model  of  33 
modes  (66  states).  A  program  flow  chart  is  shown  in  Figure  7. 

The  program  provides  decoupled  control  by  reducing  the  closed  loop 
state  equation  matrix,  given  in  form  by  Eq  (34),  to  the  lower  block  tri¬ 
angular  form.  It  is  unnecessary  to  form  the  upper  block  system  since  the 
eigenvalue  results  for  both  systems  are  identical.  In  fact,  by  renumbering 
the  controllers  in  the  spillover  elimination  portion  of  this  program,  the 
opposite  transformation  is  achieved.  For  this  reason,  only  a  main  program 
for  the  lower  block  triangular  controlled  system  is  presented.  This  is 
listed  in  Appendix  A.  The  subroutines  which  support  this  program  are  also 
in  Appendix  A.  Several  other  subroutines  are  called,  but  not  listed. 

These  are  provided  by  the  International  Mathematical  and  Statistical  Li¬ 
brary  (IMSL). 

Since  program  flexibility  is  desired,  once  the  modal  data  is  read, 
the  program  is  designed  to  make  any  number  of  runs  with  different  param¬ 
eters  for  each  run.  The  parameters  that  may  be  varied  by  the  operator  in¬ 
clude:  using  a  three  or  a  four  controller  system,  which  modes  are  assigned 
to  each  controller,  what  control  and  observer  weighting  values  to  assign 
to  each  mode,  and  what  initial  system  damping  ratio  is  applied. 

The  program  may  read  the  input  data  from  initialization  assignments 
within  the  program  or  from  a  permanent  data  file.  In  either  case,  the 
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Figure  7.  Computer  Program  Flow  Chart 


program  operates  as  if  it  were  interactive  by  prompting  for  input  and  then 
echoing  the  data  read  in.  This  makes  the  output  very  easy  to  interpret  by 
allowing  the  user  to  trace  the  computer's  progress  through  the  execution  of 
the  program. 

The  program  is  initialized  by  inputting  the  number  of  controllers  de¬ 
sired  and  then  the  number  of  modes  in  each  controller.  If  three  controllers 
are  used,  the  number  of  residual  modes  will  be  requested,  otherwise  the 
fourth  controller  system  is  run  without  residuals.  Next,  the  number  of  ac¬ 
tuators  and  sensors  are  input,  along  with  the  damping  ratio  5.  CSDL  2  was 
tested  with  twenty-one  actuators,  twenty-one  sensors,  and  a  damping  ratio 
of  0.01  for  each  mode  (15).  The  program  will  then  read  from  a  permanent 
file  the  matrix  ($TD)  of  modal  data  at  each  actuator  location,  followed  by 
the  transpose  of  the  matrix  of  modal  data  at  each  sensor  location.  In  this 
study,  these  two  matrices  are  identical  since  colocated  pairs  of  actuators 
and  sensors  are  employed.  The  sensor  modal  matrix  is  input  in  transposed 
form  so  that  the  matrix  for  a  colocated  system  may  be  copied  directly  from 
the  actuator  modal  data  matrix.  However,  these  are  left  as  separate  entries 
in  the  event  the  actuators  and  sensors  are  not  colocated.  Finally,  the  modal 
frequencies  are  read  in  from  the  data  file.  After  this  preload  of  data,  the 
desired  run  is  made  by  specifying  which  modes  are  to  be  controlled  by  each 
controller  and  which  modes  are  to  be  left  as  residuals,  along  with  the  de¬ 
sired  control  and  observer  weighting  values  for  each  mode.  Different  con¬ 
figuration  runs  may  be  made  changing  mode  assignments,  weighting  values  or 
controller  configurations. 

Program  executir  actually  begins  with  the  formation  of  the  A,  B,  C, 
and  weighting  (Q)  matrices  for  each  controller.  This  is  conveniently  done 
by  subroutines  which  read  the  required  data  for  the  modes  specified.  These 
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subroutines  allow  the  operator  to  change  the  size  of  the  controllers  as 
needed  simply  by  specifying  the  number  of  modes  to  be  placed  in  each  con¬ 
troller. 

Once  the  initial  matrices  are  formed,  the  optimal  control  gain  ma¬ 
trices,  ,  are  determined  using  a  series  of  subroutines  which  generate  a 
numerical  solution  to  the  matrix  Riccati  equation.  These  sophisticated 
routines  were  created  by  Kleinman  (7)  and  so  are  known  as  the  Kleinman 
routines.  The  matrices,  along  with  the  parameter  matrices.  A,  B,  and 
C,  are  then  combined  to  form  the  optimal  closed  loop  system  matrix.  This 
particular  program  develops  the  three  controller  system  with  residuals, 
but  does  not  include  residual  terms  in  the  four  controller  systems. 

Eigenvalue  analysis  of  the  individual  optimal  controllers  is  per¬ 
formed  next,  making  use  of  the  ISML  routine  EIGRF  which  determines  the 
eigenvalues  of  real  non-symmetric  matrices.  These  eigenvalues  are  the  de¬ 
sired  values  of  the  individual  controllers  and  are  used  to  compare  the 
suboptimal  controllers  eigenvalues  against.  The  optimal  gains  and  gener¬ 
alized  inverses  of  the  output  matrices  are  used  to  calculate  the  single 

✓s 

controller  DOFB  gains,  K.  These  suboptimal  gains  are  then  used  to  form 
closed  loop  single  controller  systems  whose  eigenvalues  are  compared  with 
the  desired  optimal  values.  The  multiple  controller  system  is  formed  with 
the  summed  gains,  <,  and  the  system  eigenvalues  are  calculated  to  evaluate 
coupling  and  use  for  comparison  to  the  transformed  system.  Next  the  de¬ 
coupling  transformation  matrices,  T  and  r,  are  obtained  and  used  to  find 
B*,  C*,  G*,  K*,  and  to  form  the  control  design  gains,  K. 

Spillover  elimination  is  the  next  step  in  the  control  algorithm. 

This  varies  with  controller  configuration.  The  3  or  4  controller  varia¬ 
tion  is  accommodated  by  the  program.  The  spillover  elimination  is  a  rather 


lengthy  portion  of  the  program.  The  modes  to  be  suppressed  are  formed  into 
non-zero  B$  and  C$^  matrices  of  the  form  given  in  Eq  (80).  The  IMSL  rou¬ 
tine  LSVDF  is  used  to  perform  a  singular  value  decomposition  on  these  ma¬ 
trices.  By  using  the  left  singular  vectors  associated  with  the  zero  sin¬ 
gular  values  of  B  and  C J ,  the  transformation  matrices  T.  and  r.  are 
formed.  The  program  then  applies  the  transformation  matrices  and  as  dis¬ 
cussed  in  Section  IV  to  create  new  matrices,  B*,  C|,  G|,  K*,  and  .  The 
new  closed  loop  systems  are  generated  and  the  eigenvalue  analysis  is  re¬ 
peated. 

This  program  demonstrates  the  flexibility  of  the  control  method  ap¬ 
plied.  The  only  change  that  needs  to  be  made  to  adapt  the  method  for  an¬ 
other  structure  is  the  basic  matrix  dimensioning  in  the  program.  Nothing 
else  has  to  be  altered,  as  long  as  the  system  model  can  be  defined  by 


x  =  Ax  +  Bu  (5) 

and 

y  =  Cx  (8) 


as  discussed  in  Section  III.  Thus,  the  ease  of  application  of  the  control 
method  described  can  be  seen.  Now,  the  performance  of  control  method  will 
be  examined  using  the  programs'  eigenvalue  analysis  results. 


VI.  Investigation 


Applying  the  decentralized  control  method  to  a  large  design  model 
is  a  complicated  task,  so  a  systematic,  building  block  approach  was  used 
to  conduct  the  study.  It  was  necessary  to  establish  the  method  of  feed¬ 
back  gain  calculation,  examine  sensor  information  configurations  (rate 
vs.  position),  and  verify  the  results  on  a  "perfect"  (no  residuals)  model. 
The  results  are  compared  against,  those  for  steady  full -state  feedback  op¬ 
timal  control.  Once  this  foundation  was  laid,  the  model  could  be  expanded 
to  the  multiple  controller  case  with  and  without  residual  modes. 

In  the  initial  phase,  a  single  controller  was  developed  to  establish 

the  method  of  feedback  gain  calculation.  The  controller  was  assigned  four 

flexible  modes  for  this  analysis.  The  four  modes  were  arbitrarily  selected 

from  the  twelve  modes  considered  critical  in  the  ACOSS  Five  report  (15), 

which  are  identified  in  Table  IV.  The  control  weighting  matrix,  Q,  was 

fixed  at  1000  I.  This  program  was  modified  to  allow  four  configurations 

of  the  output  matrix:  position  sensors  only,  C  =  [  C  .$  I  0  ]  ,  rate 

PJ 

sensors  only,  C  =  [  0  :  C  .$  ]  ,  a  mix  of  position  and  rate  sensors, 

•  VJ 

c  *  [  c®  C  $  0  0  :  0  0  C  $  C  $  ]  ,  and  both  position  and  rate 

sensors.  These  single  controller  responses  were  compared  with  those  for 
the  desired  optimal  controllers  and  to  each  other. 

Next,  the  system  was  expanded  to  the  multiple  controller  case  to  ac¬ 
tively  control  twelve  modes.  No  residuals  were  included  in  this  phase  of 
the  study.  The  twelve  modes  selected  are  the  same  twelve  identified  in 
Table  IV.  Three  and  four  controller  cases  were  examined  using  low  control 
weighting  for  the  rigid  body  modes  and  higher  (5,000  to  50,000)  weighting 


for  the  flexible  modes  as  suggested  by  Aldridge  (1).  These  controllers 
were  operated  using  either  position  sensors  only  or  rate  sensors  only. 

Mode  assignments  were  rearranged  to  study  mode  grouping,  and  control 
weighting  was  varied  to  see  the  effect  this  had  on  mode  controllability. 

Finally  the  model  was  run  using  33  modes  (4-36)  in  three  controllers 
with  residuals  and  four  controllers  without  residuals.  For  this  phase, 
only  rate  sensor  data  was  used.  In  this  expanded  model,  modal  selection 
and  grouping  were  examined  as  well  as  the  effects  of  biased  weighting  of 
lower  frequency  modes.  The  results  of  this  outline  will  be  presented 
shortly. 

One  of  the  more  difficult  tasks  in  controlling  a  large  space  struc¬ 
ture  is  the  determination  of  which  structural  modes  are  to  be  actively 
controlled  and  which  are  left  as  residuals.  Factors  which  affect  this 
selection  process  include  sensor/actuator  placement  and  alignment,  con¬ 
troller  bandwidths,  and  control  constraints  such  as  1 ine-of-sight  toler¬ 
ances. 

For  multiple  controllers,  an  additional  step  has  to  be  taken.  This 
is  the  assignment  of  the  modes  to  be  controlled  to  minimize  the  control 
effort.  Organizing  modes  into  compatible  groups  can  be  done  by  a  simple 
examination  of  the  angles  between  the  vectors  of  modal  amplitudes  (the 
angles  between  the  rows  of  the  $"*D  matrix).  Defining  these  vectors  as  , 
the  angles  may  be  found  from  the  equation  for  the  dot  product  of  two  vec¬ 
tors,  given  by 


and 


*i  *  rill^jl  coseij 


eij  =  cos 


-i 


*1  * 


*i  * 


(123) 


(124) 
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The  modes  are  then  grouped  so  no  two  orthogonal  modes  are  in  the  same  con¬ 
troller.  In  some  cases  the  grouping  of  orthogonal  modes  in  the  same  con¬ 
troller  may  be  unavoidable. 

The  selection  of  modes  based  on  angles  will  vary  with  the  number  of 
modes  in  the  design  model,  since  the  modal  vectors,  ip,  are  calculated  based 
upon  the  evaluation  model's  (all  modeled  modes)  modal  matrix.  Applying  the 
relationship  given  in  Eq  (124),  the  twelve  mode  model  relative  angles  are 
given  in  Table  V  and  for  the  36  mode  model  the  relative  angles  are  given 
in  Table  VI. 

Another  consideration  in  modal  grouping  is  unpredictable  ill  condi¬ 
tioning  of  matrices.  Aldridge  found  that  when  rigid  and  flexible  body  modes 
were  assigned  to  the  same  controller,  loss  of  symmetry  could  result  and  the 
feedback  gains  tended  to  infinity.  All  of  his  twelve  mode  controllers  re¬ 
quired  a  separate  controller  for  the  rigid  body  modes.  The  actual  groupings 
selected  for  this  study  will  be  presented  in  the  following  section. 
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VII.  Results 


As  in  most  research  efforts,  much  more  data  was  produced  than  will  be 
presented  in  this  section.  The  eigenvalue  results  that  do  appear  should  be 
considered  representative  of  the  controller  system  configurations  studied. 
Results  are  presented  in  order  of  system  magnitude,  i.e.  single  controller, 
twelve  mode  multiple  controller,  and  33  mode  multiple  controller. 

The  first  controller  studied  was  a  four  mode  single  controller.  The 
open  loop  damping  ratio  applied  was  0.01  and  a  control  weighting  matrix  of 
Q  =  1000  I_  was  used  for  modes  12,  17,  22,  and  24.  The  response  of  this 
controller  using  position  sensors  only,  rate  sensors  only,  a  mix  of  posi¬ 
tion  and  rate  sensors,  and  both  position  and  rate  sensors  (full  C  matrix) 
is  given  in  Table  VII.  Next  the  rigid  body  modes  were  placed  in  a  single 
controller  with  position  sensors  and  again  with  rate  sensors.  The  control 
weighting  matrix  was  set  at  Q  *  2  I.  The  results  for  these  controllers 
are  given  in  Table  VIII.  The  performance  of  the  single  controllers  re¬ 
vealed  that  rate  sensors  provided  the  desired  response  for  flexible  body 
modes,  but  the  response  of  the  rigid  body  modes  was  undamped  for  position 
sensor  information  and  overdamped  for  rate  sensor  information.  The  rigid 
body  modes,  however,  were  retained  in  the  12  and  33  mode  design  models,  and 
for  these  multiple  controller  runs,  only  rate  sensors  were  used. 

The  twelve  mode  model  was  examined  with  three  and  four  controllers. 
The  modes  selected  for  this  system  are  the  12  critical  modes  identified 
in  the  AC0SS  Five  report  (15),  i.e.  modes  4,  5,  6,  7,  12,  13,  17,  21,  22, 
24,  28,  30.  The  initial  mode  grouping  was  identical  to  Aldridge's  three 
controller  (1).  The  groupings  are: 
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Controller  1:  4,  5,  6 

Controller  2:  7,  13,  17,  21,  30 

Controller  3:  12,  22,  24,  23 

The  groupings  in  these  controllers  appear  undesirable  since  modes  pairs 
4-5,  4-6,  7-13,  7-17,  7-30,  13-21,  17-21,  21-30,  12-28,  22-28,  and  24-28 
are  essentially  orthogonal  (see  Table  V).  The  control  weighting  matrix 
for  the  flexible  modes  was  set  at  Q  =  5000  I_  and  remained  at  Q  =  2 
for  the  rigid  body  modes.  The  eigenvalues  for  the  closed  loop  system  ap¬ 
pear  in  Table  IX.  While  the  untransformed  system  shows  improved  damping 
for  six  of  the  nine  flexible  modes,  the  controller  coupling  has  driven 
mode  13  (3.74  rad  s-1)  unstable.  After  the  transformation,  the  flexible 
modes  remain  stable  with  imp.  oved  closed  loop  damping.  Modes  28  and  30 
remain  unchanged  though.  These  two  nearly  orthogonal  modes  were  subse¬ 
quently  given  a  higher  control  weighting  and  reassigned  to  other  mode 
groupings: 


Controller  1: 

4,  5,  6 

Controller  2: 

7,  13,  17,  21 

Controller  3: 

12,  22,  24,  28,  30 

Controller  1: 

4,  5,  6 

Controller  2: 

7,  13,  17,  21 

Controller  3: 

12,  22,  24 

Controller  4: 

28,  30 

However,  the  closed  loop  damping  of  modes  28  and  30  remained  unchanged  from 
the  open  loop  damping  of  .010.  These  two  modes  appear  weakly  controllable 
with  the  given  actuator  configuration.  Examination  of  the  modal  matrix 
(Appendix  B)  shows  the  elements  for  modes  28  and  30  are  zero  or  "small" 


compared  to  the  elements  for  the  other  modes,  which  supports  the  condition 
of  weakly  control! able/observable.  Also  of  note  is  that  orthogonal  mode 
pairs  like  7-13  not  only  remained  controllable,  but  both  show  nice  improve¬ 
ments  of  closed  loop  damping.  The  optimal  control  performance  before  and 
after  transformation  is  also  of  interest  and  is  compared  for  this  system 
in  Table  X. 

With  the  multiple  controller  program  operating,  the  next  step  was  to 
apply  it  to  the  33  mode  design  model.  Aldridge's  results  showed  ill  con¬ 
ditioning  of  matrices  when  the  rigid  body  modes  were  grouped  with  flexible 
modes.  So  for  the  first  run  the  rigid  modes  were  grouped  in  a  separate 
controller  and  the  system  modes  were  grouped  as: 

Controller  1:  4,  5,  6 

Controller  2:  7,  8,  9,  10,  11,  12,  13,  14,  15,  16 

Controller  3:  17,  18,  19,  20,  21,  22,  23,  24,  25,  26 

Controller  4:  27,  28,  29,  30,  31,  32,  33,  34,  35,  36 

The  eigenvalue  results  for  this  system  are  given  in  Table  XI.  Here,  as 

with  the  12  mode  case,  the  untransformed  system  shows  coupling  that  de¬ 

stabilizes  modes  8,  18,  20,  21,  and  30,  with  mode  8  driven  unstable.  The 
transformed  system  remains  completely  stable  with  only  one  mode  (11)  suf¬ 
fering  some  damping  loss.  Most  modes  closed  loop  damping  improves.  Es¬ 
pecially  note  controller  3  (modes  17-26)  whose  closed  loop  damping  for  all 
inodes  has  at  least  doubled.  Five  of  the  higher  frequency  modes  in  con¬ 
troller  4  improved,  including  mode  30  which  was  uncontrollable  in  the  12 
mode  model.  The  modes  whose  closed  loop  damping  remained  equal  to  the 
passive  .010  value  may  again  be  weakly  observable/controllable  with  the 
given  actuator/sensor  configuration. 
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At  this  point  there  were  two  system  conditions  to  investigate.  First 
was  the  ill  conditioning  of  matrices  that  was  seen  when  grouping  rigid  and 
flexible  modes  in  one  controller.  And  second,  was  the  initial  33  mode 
grouping  denying  better  control  of  the  higher  frequencies.  To  check  for 
ill  conditioning,  the  rigid  body  modes  were  mixed  with  flexible  modes  in 
four  different  configurations  as  follows: 

Grouping  1:  4,  5,  6,  7,  9,  10,  13 
Grouping  2:  4,  5,  6,  15,  16,  17,  18,  19 

Grouping  3:  4,  5,  6,  23,  24,  25,  26,  27 

Grouping  4:  4,  5,  6,  31,  32,  33,  34,  35,  36 

Mixing  with  groupings  1,  2,  and  3  resulted  in  completely  stable  decoupled 

controllers,  but  when  grouped  with  the  higher  frequency  modes  (Grouping  4), 
the  controller  gains  became  large  and  the  eigenvalue  analysis  was  stifled 
with  infinite  values.  To  see  if  this  was  not  just  a  problem  of  mixing  rigid 
modes  with  the  higher  frequencies,  the  lower  frequency  modes  were  mixed 
with  the  higher  ones  as  well: 

Controller  1:  4,  5,  6,  15,  16,  17,  18,  19 

Controller  2:  7,  8,  9,  23,  24,  33,  34,  35 

Controller  3:  10,  11,  12,  25,  26,  31,  32,  36 

Controller  4:  13,  14,  20,  21,  22,  27,  28,  29,  30 

This  configuration  was  "nice"  and  eigenvalue  results  after  transformation 
are  given  in  Table  XII.  Mode  28  remained  the  only  mode  whose  closed  loop 
damping  could  not  be  improved  via  modal  grouping.  Unfortunately  this  mode 
is  considered  critical  to  LOS  performance. 

Finally  the  33  mode  model  was  run  in  the  three  controller  configura¬ 
tion  with  residuals.  Most  of  the  configurations  that  contained  more  than 
five  residuals  drove  one  or  more  of  the  transformed  system  modes  unstable. 


However,  two  modal  groupings  that  were  successful,  i.e.  remained  stable 
after  transformation,  were  one  with  three  residual  modes  and  the  other 
with  seven  residuals.  The  first  system  retained  weighting  matrices  of 

Qrigid  =  2  I  and  ^flexible  =  5000  -  and  was  grouped  as: 

Controller  1:  4,  5,  6,  9,  10,  11,  28,  29 

Controller  2:  8,  12,  14,  15,  16,  17,  18,  19,  30,  31 
Controller  3:  20,  21,  22,  23,  24,  25,  26,  27,  32,  33 
Residuals:  34,  35,  36 

The  system  eigenvalue  results  for  this  configuration  are  given  in  Table 
XIII.  For  this  configuration,  mode  28  grudgingly  improved  and  most  of  the 
actively  controlled  modes  show  improved  closed  loop  damping. 

The  last  controller  system  to  examine  has  the  modal  groupings: 
Controller  1:  4,  5,  6,  9,  10,  11,  13 

Controller  2:  8,  12,  14,  15,  16,  17,  18,  19,  20 

Controller  3:  21,  22,  23,  24,  26,  28,  29,  30,  31 

Residuals:  25,  27,  32,  33,  34,  35,  36 
The  control  weighting  matrix  on  modes  7-20  was  raised  to  Q  =  50,000  I_ 
and  lowered  on  modes  21-36  to  create  high  control  of  mode  frequencies  be¬ 
low  6  rad  s"1,  and  a  "dead  zone"  effect  for  the  frequencies  above  that 
cutoff.  The  eigenvalue  results  for  this  system  are  given  in  Table  XIV. 
Again  the  closed  loop  damping  is  remarkably  improved  for  modes  20  and  be¬ 
low  and  only  the  residual  mode  35  loses  some  damping. 

In  general,  the  untransformed  system  shows  the  destabilizing  effect 
of  spillover.  After  transformation,  though,  system  performance  still  is 
highly  dependent  upon  modal  grouping  and  control  weighting.  The  condi¬ 
tion  and  degree  of  controllability/observability  is  absolutely  dependent 
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on  the  actuator/sensor  configuration  and  modal  matrix  values.  This  is 
easily  seen  in  following  the  response  of  mode  28  through  it  various  con¬ 
troller  assignments. 


TABLE  VII 


Eigenvalue  Analysis  -  Single  Controller 
Modes:  12,  17,  22,  24 


OPTIMAL  FULL  STATE  FEEDBACK  ? 


-.6978118  +  11. 1193805i  .063 
-.1950110+  7. 2780866 i  .027 
-.0613474+  5. 1212778i  .012 
-.1814128+  3.4979155i  .052 

DOFB  POSITION  SENSORS 

-.111409  +  11. 1441 i  .010 
-.072806  +  7. 2804 i  .010 
-.051216  +  5. 12 14i  .010 
-.035193  +  3.5013i  .010 

DOFB  RATE  SENSORS 

-.697751  +  11. 115567 i  .063 
-.194998  +  7.277893i  .027 
-.061347  +  5. 121270i  .012 
-.181487  +  3.498374i  .052 

POSITION  (12.17)  AND  RATE  (22,24) 

-.697816  +  11 . 117928i  .063 
-.194983  +  7.278735i  .027 
-.051216  +  5. 1213891  .010 
-.035056  +  3.501293i  .010 

POSITION  AND  RATE  SENSORS 

-.425969  +  11 . 160118i  .038 
-.134633  +  7.287732i  .018 
-.056303  +  5. 122321i  .011 
-.107221  +  3 . 52 12 15i  .030 
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TABLE  VIII 

Eigenvalue  Analysis  -  Single  Controller 
Rigid  Body  Modes:  4,  5,  6 


OPTIMAL 

-0.16093  +  . 15691 i 
-0.14062  7  . 13792i 
-0.11454  7  . 11307 i 


POSITION  SENSORS 

0.0  +  . 2248i 
0.0  +  . 1970i 
0.0  +  . 1609i 


RATE  SENSORS 


-.55906  +  O.Oi  \ 
-.00018  +  O.Oi  J 

»1 

-.17413  +  O.Oi  \ 
-.00057  +  O.Oi  i 

»1 

-.04497  +  O.Oi  \ 
-.00222  +  O.Oi  f 

»1 

TABLE  IX 


Eigenvalue  Analysis  -  3  Controllers 

Q  .  . .  *  2  Q..  ...  =  5000 

yrigid  ^flexible 


Mode  Assignments 

Controller  1:  4,  5,  6  Controller  3:  12,  22,  24,  28 

Controller  2:  7,  13,  17,  21,  30  Residuals:  None 


Overall  System  Eigenvalues 


Mode 

Before  Transformation 

After  Transformation 

1 

30 

-0.2505  +  25.0483i 

.010 

-0.2505  +  25. 0483 i 

.010 

28 

-0.2169  +  21.69031 

.010 

-0.2169  +  21.69031 

.010 

24 

-2.7170  +  10. 52901 

.250 

-0.7316  + 

11.11061 

.066 

22 

-0.3863  + 

7.1399i 

.054 

-0.2549  + 

7.26531 

.035 

21 

-0.2484  + 

6.10621 

.041 

-0.2987  + 

6.09881 

.049 

17 

-0.0679  + 

5.12491 

.013 

-0.0853  + 

5. 11911 

.017 

13 

+0.3818  + 

4.1764i 

unstable 

-0.9443  + 

3.62711 

.252 

12 

-0. 1040  + 

3.3917i 

.031 

-0.2371  + 

3.50111 

.068 

7 

-0.0191  + 

0.6893i 

.028 

-0.0427  + 

3.71481 

.060 

6 

-7.1082  + 

0.01  \ 

»1 

-0.3219  + 

O.Oi  \ 

»1 

0.1  E-8  + 

0.01  f 

0.0002  + 

O.Oi  f 

5 

-2.1548  + 

O.Oi  \ 

»1 

-0.2535  + 

O.Oi  \ 

»1 

0.2  E-7  + 

O.Oi  f 

0.0002  + 

O.Oi  ) 

4 

-0.4526  + 

O.Oi  l 

»1 

0.1414  + 

O.Oi  \ 

»1 

0.0001  + 

O.Oi  f 

0.0023  + 

O.Oi  f 

61 


*i  i 


TABLE  X 


& 

V& 

A 

A 


Eigenvalue  Analysis  -  3  Controllers  (Optimal) 
Qr1gid  =  2  ^flexible  =  5000 


Mode  Assignments 

Controller  1:  4,  5,  6  Controller  3:  12,  22,  24,  28 

Controller  2:  7,  13,  17,  21,  30  Residuals:  None 


Controller  Eigenvalues 

Before  Transformation  (BG)  £  After  Transformation  (B*G*)  £ 


Controller  1 


-0.16093  +  0.1569H 

.716 

-0.16093  +  0.1569H 

.716 

-0.14062  +  0.137921 

.714 

-0.14062  +  0. 137921 

.714 

-0.11454  +0.113071 

.712 

-0.11454  +  0.113071 

.712 

Controller  2 

-0.25050  +  25.04834i 

.010 

-0.25050  +  25. 04834 i 

.010 

-0.34700  +  6.098071 

.057 

-0.29920+  6.100411 

.049 

-0.08782  +  5.120711 

.017 

-0.08661  +  5.120821 

.017 

-1.12614+  3.596201 

.299 

-0.94251  +  3. 641741 

.251 

-0.59706  +  0.755481 

.620 

-0.04268  +  0.716441 

.059 

Controller  3 

-0.21692  +  21 . 69033i 

.010 

-0.21691  +  2 1.69033 i 

.010 

-1.54350  +  11.034651 

.139 

-0.73205  +  11. 11716i 

.066 

-0.40996  +  7. 26962 i 

.056 

-0.25523  +  7.276241 

.035 

-0.39925  +  3.482631 

.114 

-0.23632  +  3.49517 1 

.067 

«V«\ 
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TABLE  XI 


Eigenvalue  Analysis  -  4  Controllers 
Mode  Assignments 


Controller  1:  4-6 
Controller  2:  7-16 


Controller  3:  17-26 

Controller  4:  27-36 


Before  Transformation 


After  Transformation 


-0.5563 
-0.5311 
-0.5066 
-0.4952 
-2.1740 
-0.3027 
-0.2300 
-0.2486 
-0.2170 
-0.2274 
-0.1372 
-0.2164 
-0.1203 
-0. 1085 
-0.1823 
-0.0529 
-0.0463 
-0.0535 
-0.0093 
-0.0987 
-0.0511 
-0.0409 
-0.0643 
-0.0855 
-0.1230 
-0.1786 
-0.0247 
-0.0100 
+0.0242 
-0.4161 


+  55. 6271 1 
+  53.1101i 
+  50. 6708 1 
+  41 . 2112i 
+  26. 8590i 
+  25. 5128i 
+  25. 0623i 
+  24. 8632 i 
+  21. 6915i 
+  21.57691 
+  14. 1619i 
+  14. 12631 
+  11.13741 
+  9.74191 
+  7.2325i 
+  6.10241 
+  5.67001 
+  5.1690i 
+  5 . 1363i 
+  5. 1204 i 
+  4.00701 
+  3. 9980 i 
+  3.8315i 
+  3.72211 
+  3. 4744 i 
+  2.81241 
+  1.00291 
+  0.9235i 
+  0.72611 
+  0.66631 


-0.9692  +  0.01 
-0.0002  +  O.Oi 

-0.1755  +  O.Oi 
-0.0006  +  O.Oi 

-0.0502  +  O.Oi 
-0.0022  +  O.Oi 


TABLE  XII 


Controller  1: 
Controller  2: 


Eigenvalue  Analysis  -  4  Controllers 


Mode  Assignments 

4-6,  15-19  Controller  3:  10-12,25,26,31,32,36 

7-9,23,23,33,34,35  Controller  4:  13,14,20-22,27-30 


After  Transformation 


Controller  1 

£ 

-0.0590  + 

5.172H 

.011 

-0.0513  + 

5. 1248 i 

.010 

-0.4740  + 

4.9394i 

.096 

-0.1130  + 

4.0306i 

.028 

-0.0406  + 

3.9992i 

.010 

Rigid  Body 

»1 

Controller  2 

£ 

-0.6904  + 

53. 1049i 

.013 

-0.5774  +  50. 6682 i 

.011 

-0.4449  +  41 . 23141 

.011 

-0.3913  + 

1 1 . 08  Hi 

.035 

-0.3550  + 

9.7809i 

.036 

-0.3861  + 

0.3432i 

.416 

-0.4367  + 

0.61731 

.578 

-0.7272  + 

0.5072i 

.820 

Controller  3 

£ 

-0.5927  +  55.62531 

.011 

-0.4866  +  27. 2171i 

.018 

-0.5558  +  25. 4804 i 

.022 

-0.1890  + 

14. 1351i 

.013 

-0.9301  + 

14.09281 

.066 

-0.1337  + 

3.4590i 

.039 

-0.3617  + 

2. 8683 i 

.125 

-0.5185  + 

0.9766i 

.469 

Controller  4 

£ 

-0.2513  +  25.0480i 

.010 

-0.3108  +  24.8594i 

.013 

-0.2210  +  21 . 6878i 

.010 

-0.2834  +  21.57321 

.013 

-0.1345  + 

7.27051 

.018 

-0.0616  + 

6. 1070i 

.010 

-0.1522  + 

5. 7430 i 

.026 

-0.0564  + 

3.80721 

.015 

-0.2251  + 

3. 8044 i 

.059 

TABLE  XIII 


Controller  1 
Controller  2 


Eigenvalue  Analysis  -  3  Controller 


Mode  Assignments 

:  4-6,  9-11,  28,  29  Controller  3:  20-27,  32, 

:  8,  12,  14-19,  30,  31  Residuals:  34-36 


Mode 

After  Transformation 

£ 

36 

-0.5562  +  55.6270i 

.010 

35 

-0.5313  + 

53 . 1227 i 

.010 

34 

-0.5067  + 

50. 6687 i 

.010 

33 

-0.4124  +  41 . 2324i 

.010 

32 

-0.2730  +  27 . 2535i 

.010 

31 

-0.6159  + 

25. 4194i 

.024 

30 

-0.6987  +  25 . 0474i 

.028 

29 

-0.2496  +  24. 8631 i 

.010 

28 

-0.2484  +  21. 6868 i 

.011 

27 

-0.2164  +  21 . 5727 i 

.010 

26 

-0.2615  +  14. 1483i 

.018 

25 

-0.1416  +  14. 1373i 

.010 

24 

-0.2567  + 

11 . 1296i 

.023 

23 

-0.1845  + 

9.7514i 

.019 

22 

-0.0785  + 

7. 2801 i 

.011 

21 

-0.0864  + 

6. 1025i 

.014 

20 

-0.1328  + 

5.7604i 

.023 

19 

-0.6081  + 

5 . 0370i 

.120 

18 

-0.4594  + 

5 . 0346i 

.091 

17 

-0.2799  + 

5 . 0143i 

.056 

16 

-0.7299  + 

4. 0903 i 

.176 

15 

-0.2282  + 

3. 9903 i 

.057 

14 

-0.7534  + 

3.7874i 

.195 

13 

-0.0503  + 

3.7379i 

.013 

12 

-0.2638  + 

3. 5546 i 

.074 

11 

-0.0310  + 

2. 8596 i 

.011 

10 

-0.0149  + 

1.0883i 

.014 

9 

-0.0364  + 

0.8697i 

.042 

8 

-0.0183  + 

0.6945i 

.026 

7 

-0.7512  + 

0.5480i 

.808 

4-5-6 

»1 

TABLE  XIV 

Eigenvalue  Analysis  -  3  Controllers 


Controller 

Controller 


Mode  Assignments 

1:  4-6,  9-11,  13  Controller  3:  21-24,  26,  28-31 
2:  8,  12,  14-20  Residuals:  25,  27,  32-36 


Mode 

After  Transformation 

36 

-0.5488  +  55. 6505i 

.010 

35 

-0.3455  +  53. 1892 i 

.006 

34 

-0.5839  +  50.6627i 

.012 

33 

-0.5755  +  41.0474i 

.014 

32 

-2.0730  +  26 . 9359i 

.077 

31 

-0.2586  +  25. 4573i 

.010 

30 

-0.2524  +  25. 0479i 

.010 

29 

-0.2507  +  24.8623i 

.010 

28 

-0.2172  +  21. 6903 i 

.010 

27 

-0.4040  +  21.5123i 

.019 

26 

-0.1450  +  14.16251 

.010 

25 

-0.2526  +  14. 1 1391 

.018 

24 

-0.1137  +  11.14041 

.010 

23 

-0.0989  +  9. 7452 i 

.010 

22 

-0.0740  +  7. 2802 i 

.010 

21 

-0.0643  +  6.1075i 

.011 

20 

-0.3642  +  4.97811 

.073 

19 

-0.3872  +  4. 8587 i 

.079 

18 

-1.2217  +  4. 8231 i 

.246 

17 

-1.9930  +  4.2648i 

.423 

16 

-2.2082  +  3.9841i 

.485 

15 

-0.9154  +  3.92371 

.227 

14 

-0.4845  +  3.7452i 

.128 

13 

-0.1403  +  3.58731 

.039 

12 

-2.6692  +  3.32071 

.627 

11 

-3.0419  +  2.63711 

.756 

10 

-0.2659  +  2.38181 

.111 

9 

-0.0274  +  0.85651 

.032 

8 

-0.1131  +  0.25241 

.409 

7 

-0.1027  +  0.03861 

.936 

-5-6 

»1 
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VIII.  Conclusions 


This  investigation  demonstrated  the  feasibility  and  merits  of  the 
decentralized  direct  output  feedback  control  method.  This  technique  can 
be  used  to  design  feedback  gains  for  multiple  controllers  that  decouple 
the  individual  controllers  by  eliminating  controller  to  controller  obser¬ 
vation  (or  control)  spillover.  While  system  stability  was  maintained  for 
design  models  that  did  not  include  residual  modes,  system  stability  can¬ 
not  be  guaranteed  for  design  models  that  include  residual  modes.  Modal 
grouping  played  a  significant  part  in  the  system  stability  achieved. 

Suitable  modal  grouping  may  have  to  be  accomplished  by  "trial  and 
error"  unless  a  more  technical  approach  can  be  devised.  The  use  of  angles 
between  modal  amplitude  or  rate  vectors  is  a  convenient  method  for  initial 
grouping  of  modes.  However,  orthogonal  modal  pairs  may  be  satisfactorily 
controlled  even  though  they  are  in  the  same  controller.  The  rank  of  the 
B  and  C  matrices  should  be  examined  closely.  If  they  are  not  of  full  rank, 
the  modal  groups  may  not  be  fully  compatible  as  indicated  by  zero  entries 
in  the  non-zero  singular  values  of  the  singular  value  decomposition.  Since 
the  B  and  C  matrices  are  a  function  of  actuator  and  sensor  placement,  this 
model  evaluation  technique  is  a  valuable  design  tool  for  examining  perfor¬ 
mance  of  actuator  and  sensor  configurations.  The  matrix  can  be  checked 
for  determining  mode  observability  and  controllability.  This  step  should 
be  considered  part  of  the  control -configured  system  approach. 

A  major  shortcoming  of  direct  output  feedback  using  rate  sensors  was 
the  loss  of  control  of  the  rigid  body  modes.  This  is  certainly  not  a  fatal 
flaw  for  this  control  method.  The  rigid  body  modes  may  be  controlled  by  a 
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separate  optimal  control  system  or  additional  states  (e.g.,  rigid  body  an¬ 
gles)  may  be  added  to  the  state  vector. 

Overall  closed  loop  damping  was  improved  for  actively  controlled 
modes  and  system  stability  maintained  for  both  the  12-  and  33-mode  design 
models.  This  study  has  demonstrated  a  control  technique  that  eliminated 
the  requirement  for  the  observer  states  of  modern  control  and  allows  either 
control  of  twice  as  many  modes  or  dramatically  reduces  the  computation  bur¬ 
den  to  control  the  same  number  of  modes. 


IX.  Recommendations 


While  this  study  has  demonstrated  the  successful  application  of  the 
decentralized  DOFB  method,  it  is  not  an  end  in  itself.  Before  a  control 
method  is  applied  to  a  system,  a  control -configuration  approach  should  be 
developed  to  establish  desirable  mode  selection  and  grouping.  This  should 
involve  examination  of  techniques  to  orient  sensors  and  actuators  specifi¬ 
cally  for  the  critical  modes.  This  approach  should  include  examination  of 
the  changes  in  the  <i>^D  matrices  with  variations  of  the  number  of  modes  in 
the  design  model  and  reconfiguration  of  sensors/actuators.  It  may  be  pos¬ 
sible  to  develop  a  method  to  establish  optimal  modal  grouping. 

System  robustness  is  a  desirable  characteristic  for  control  systems. 
Just  how  sensitive  this  control  method  is  to  parameter  variations  and  com¬ 
ponent  failure  should  be  established.  Also,  analysis  to  predict  ill  con¬ 
ditioning  of  matrices  would  certainly  be  useful  to  a  control  designer. 

The  systems  studied  have  provided  a  take  off  point  for  calculations 
of  the  line  of  sight  (LOS)  error  for  both  forced  and  unforced  responses. 

The  systems  should  also  be  modified  to  provide  proper  control  of  the  rigid 
body  modes  if  they  are  included  in  the  design  model.  This  may  be  easily 
achieved  by  addition  of  body  angles  to  the  state  vector. 

One  further  item  of  interest  is  the  computer  burden  comparisons  for 
modern  (cost  optimal)  and  DOFB  calculations.  Both  systems  should  be  written 
using  efficient  programming  techniques. 


nnnrit' 


PROGRAM  D0FB33( INPUT , OUTPUT,  TAPES, TAPES) 

C 

PEN:  D0FB335 

THIS  PROGRAM  GENERATES  A  LOWER  TRIANGULAR  TRANSFORMATION 
THE  THREE  OR  FOUR  CONTROLLER  SOLUTION  MAY  INCLUDE  RESIDUALS. 

REAL  Al(21,21), A2(21, 21), A3(21, 21), A4(21,21) 

REAL  B1 ( 21 , 21 ) , B2( 21 , 21 ) , B3(21, 21 ) , B4( 21, 21 ) 

REAL  BSTAR2 (21,21), BSTAR3 (21,21), BSTAR4 (21,21) 

REAL  CSTAR1 (21,21), CSTAR2 (21,21), CSTAR3 ( 21,21) 

REAL  Cl(21, 21), C2( 21, 21),C3(21,21), C4(21, 21 ) 

REAL  CTCC1 (21,21), CTCC2 (21,21), CTCC3 (21,21) 

REAL  GT1 (21,21), GT2 (21,21), GT3 (21,21) 

REAL  SI ( 24) , KSTAR4( 21 , 21 ) , GSTAR4( 21 , 21 ) 

REAL  UT ( 21 , 21 ) , VP(21, 21 ) , QPLUS( 21, 21 ) 

REAL  QUT ( 21 , 21 ) , CPLUS( 21 , 21 ) , KHAT4 (21,21) 

REAL  KHAT1 (21 , 21 ) , KHAT2( 21 , 21 ) , KHAT3(21, 21 ) 

REAL  KSTAR1 (21,21), KSTAR2 (21,21), KSTAR3 (21,21) 

REAL  SAT (21, 21 ) , SAT2(21 , 21 ) , SAT3( 21 , 21 ) , SAT4(21 , 21 ) 

REAL  AKC(21 , 21 ) , ACT (21 , 21 ) , BCG( 21 , 21 ) , KCC ( 21 , 21 ) 

REAL  P(21,21),S(21,21) 

REAL  ABK1(21, 21), ABK2(21,21), ABK3(21, 21) 

REAL  GSTAR1 (21,21), GSTAR2 (21,21), GSTAR3 (21,21) 

REAL  TKG1 (21, 21 ) , TKG2(21, 21 ) , TKG3(21 , 21 ) 

REAL  QA1 (21,21), QA2 (21,21), QA3( 21 , 21 ) , QA4( 21 , 21 ) 

REAL  ABG1 (21,21), ABG2 ( 21 , 21 ) , ABG3 ( 21 , 21 ) , ABG4 ( 21 , 21 ) 

REAL  ABG5( 21 , 21 ) , ABG6( 21 , 21 ) , A8G7( 21 , 21 ) , ABG8( 21 , 21 ) 

REAL  GAIN1(21, 21 ), GAIN2(21, 21), GAIN3(21, 21), GAIN4(21, 21) 
REAL  KT1 (21, 21 ) , KT2(21, 21 ) , KT3(21, 21 ) , KT4(21, 21 ) 

REAL  KOB1 (21,21), K0B2 (21,21), K0B3 ( 21 , 21 ) , KOB4 ( 21 , 21 ) 

REAL  GAMMA1(21, 21), GAMMA2(21,21), GAMMA3(21,21) 

REAL  T2(21, 21 ) , T3(21, 21 ) , T4(21, 21 ) 

REAL  TRT(21,21),TEN(21,21),CT(21,21), V(21,21) 

REAL  RK ( 21 , 21 ) , RK1 (21,21), RK2 (21,21), RK3 (21,21) 

REAL  RG2 (21,21), RG3 (21,21),  , 21 ) 

REAL  MAJM( 66, 66 ) , D( 33 ) , XO ( 66 ) ,  W( 33 ) , TOL , DT 
REAL  ZETA, AA( 33) , BB(33 ) , SING( 21 ) , XTR( 21 , 21 ) , XI (66) 

REAL  EAT ( 1 , 1 ) , EAT2 ( 1 , 1 ) , WORK ( 66 , 66 ) , STOR ( 21 , 2 1 ) 

REAL  PHIA(33, 21) , PHIS(33, 21) , MODE(2, 21), INIT(4, 21) 
INTEGER  N,  N2, NCI, NC2, NC3, NC12, NC22, NC32,  NR, NR2 
INTEGER  IC1 (21 ) ,  IC2(21 ) ,  IC3(21 ) , I, J, K, L, M, KK, LL, MM 
INTEGER  DEC, Q, NACT, NSEN, IR(21 ) , IER, SKIP, NCOL, NCOL1 
INTEGER  NDA, NDIM,  NDA1,  NDIM1,  ZZ,  E2,  E3,  E4,  P1,P2,  P3, P5 
COMPLEX  Z(24), Wl(12) 

COMMON/MAINA/NDA, NDA1 
COMMON/MAINB/NCOL, NCOL1 
COMMON/MAIN1 / NDIM, NDIM1 , TEN, X( 3364) 

COMMON/MAIN2/STOR 
COMMON/MAIN3/XTR 
COMMON/SAVE/T( 100) , TS( 100) 

COMMON/ INOUT/KOUT, TAPE 
COMMON/NUM/IC1, IC2, IC3,  IR,  NCI,  NC2,  NC3,  NR 
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c 

c 

C  INITIALIZATIONS 

C 

C 

NDIM  =  33 
N0IM1  =  34 
NCOL  =  21 
NC0L1  =  22 
NDA  =  66 
NDA1  =  67 
KOUT  =  6 
TAPE  =  9 
Q  =  0 
IER  =  0 
ZZ  =  0 
C 
C 


PRINT' (////)' 

PRINT*, •  ****************** ****** ********* »**< 

PPTMT*  * 

PRINT*] '  *****  DIRECT  OUTPUT  FEEDBACK 

PRINT*, '  ***** 

PRINT*, '  *****  34L0WER-RESIDUAL 

PRINT*, '  *****  BLOCK 

PRINT*, '  *****  C  S  D  L  II 

PRINT*, '  ***** 

PRINT*,  1  ■»>»*»*»!>■  »m»iimiii»»»ii»»innini«iiiii<»mi' 

PRINT' (//) ' 

PRINT*, '  THIS  PROGRAM  GENERATES  A  SOLUTION', 

*  '  USING  A  LOWER  TRIANGULAR  TRANSFORMATION 


PRINT' (////)' 


INITIAL  SELECTION  FOR  THREE  OR  FOUR  CONTROLLERS 


PRINT*, '  FOR  A  THREE  CONTROLLER  RUN,  ENTER  3,  OR,  ' 
PRINT*, '  FOR  A  FOUR  CONTROLLER  RUN,  ENTER  4  > ' 

R£AD(8, *)  DEC 

DEC  DEFAULT  SWITCH 

IF  (DEC.NE.4)  DEC  =  3 
PRINT*,'  ' 

PRINT*, '  THIS  IS  A  '.DEC, '  CONTROLLER  RUN  ' 

C 

C 
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C  PHI  MATRICES  ANO  CONTROLLER  ENTRIES 

C 

C 

PRINT* (///)* 

IF  (DEC.EQ.3)  THEN 

PRINT*,'  ENTER  NCI, NC2.NC3, NR, NACT, NSEN, ZETA  >■ 
ELSE 

PRINT*,'  ENTER  NCI, NC2,NC3,NC4, NACT, NSEN, ZETA  >' 
ENDIF 

READ(8,*)  NCI, NC2,NC3,  NR, NACT, NSEN, ZETA 
PRINT* , NCI , NC2 , NC3 , NR , NACT , NSEN, ZETA 
PRINT*, '  ' 

PRINT*,'  ENTER  THE  '.NACT,'  ELEMENTS  FOR  EACH  PHIA 
PRINT*,'  ' 

N  =  NCI  +  NC2  +  NC3  +  NR 
DO  1  1=1, N 

PRINT*, 'ENTER  PHIA  ',1, '  >' 

READ(8,*)  (PHIA(I.O), 0=1, NACT) 

PRINT*,'  ' , (PHIA(I, J), J=l, NACT) 

1  CONTINUE 
PRINT' (//)' 

PRINT*,'  ENTER  THE  ' , NSEN, '  ELEMENTS  FOR  EACH  PHIS 
PRINT*,'  ' 

DO  2  1=1, N 

PRINT*, 'ENTER  PHIS  ',1,'  >' 

READ(8,*)  (PHIS(I,0), J=1,NSEN) 

PRINT*,'  ',  (PHIS(I, J), J=1,NSEN) 

2  CONTINUE 
PRINT' (//) ' 

C 

C 

C  OMEGAS 

C 

C 

PRINT*,  '  ENTER  THE  VALUE  FOR  EACH  OMEGA  ' 

PRINT*, '  ' 

DO  3  1=1, N 

PRINT*, 'ENTER  OMEGA  ',1,'  >' 

READ{8,*)  W(I) 

PRINT*,'  ' , W(I) 

D(I)  =  -2.  *  ZETA  *  W(I) 

3  CONTINUE 
C 

C 

20  CONTINUE 
C 
C 

C  SECONDARY  SELECTION  FOR  THREE  OR  FOUR  CONTROLLERS,  TO 
C  BE  USED  FOR  RUNS  AFTER  THE  FIRST  JOB 
C 
C 

IF  (Q.EQ.2)  THEN 

PRINT*, '  FOR  A  THREE  CONTROLLER  RUN,  ENTER  3,  OR,  ' 
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PRINT*, 1  FOR  A  FOUR  CONTROLLER  RUN,  ENTER  4 
READ(8, *)  DEC 
C 

C  DEC  DEFAULT  SWITCH 
C 

IF  (DEC.NE.4)  DEC  =  3 
PRINT*,'  ’ 

PRINT*, '  THIS  IS  A  '.DEC, '  CONTROLLER  RUN  ' 
PRINT' (//) ' 

C 

C 

IF  (DEC.EQ.3)  THEN 

PRINT*,'  ENTER  THE  VALUES  OF  NC1,NC2,NC3,NR  >' 
ELSE 

PRINT*, '  ENTER  THE  VALUES  OF  NCI, NC2, NC3, NC4  >' 
ENDIF 

READ(8,*)  NC1,NC2,NC3,NR 
PRINT*, NCI, NC2,NC3,  NR 
PRINT' (//) ' 

ENDIF 


(t 


C 

C 

C 

c 

c 


c 

c 


MODE  ASSIGNMENT  TO  CONTROLLERS 


PRINT' (///)' 

PRINT*,  '  THE  FOLLOWING  MODES  ARE  ENTERED  ACCORDING  TO  THE 
PRINT*, '  ORDER  IN  WHICH  THEY  ARE  ENTERED  IN  THE  DATA  FILE 
PRINT*, '  AND  NOT  ACCORDING  TO  THEIR  ACTUAL  MODE  NUMBER.  ' 
PRINT' (//)' 

PRINT*,'  ENTER  THE  ',NC1,'  CONTROLLER  1  MODES  >' 

READ(8,*>  (IC1(I), I=1,NC1) 

PRINT*,*  ' , (IC1(I), 1=1, NCI) 

PRINT*, '  ' 

PRINT*, '  ENTER  THE  ',NC2, '  CONTROLLER  2  MODES  >' 

READ(8,*)  (IC2(I), IS1,NC2) 

PRINT*,'  1 , (IC2(I), I=1,NC2) 

PRINT*,'  ' 

PRINT*,'  ENTER  THE  \NC3, '  CONTROLLER  3  MODES  >’ 

READ(8,*)  (IC3(I), 1=1, NC3) 

PRINT*,'  ', (IC3(I),I=1,NC3) 

PRINT*, '  ' 

IF  (DEC.EQ.3)  THEN 

PRINT*, '  ENTER  THE  ' , NR, '  RESIDUAL  MODES  >' 

ELSE 

PRINT*, '  ENTER  THE  ' , NR, '  CONTROLLER  4  MODES  >' 

ENDIF 

READ(8,*)  ( IR( I) , 1=1, NR) 

PRINT*,'  ' , (IR(I), 1=1, NR) 

PRINT*,'  ' 


NC12  =  2  *  NCI 
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NC22  =  2  *  NC2 

NC32  =  2  *  NC3 

N2  =  2  *  N 

NR2  =  2  *  NR 

M  =  NC12  +  NC22  +  NC32  +  NR2 
NDA  =  M 
N0A1  =  M  +  1 
NDIM  =  N 
NDIM1  =  N  +  1 
100  CONTINUE 

PRINT*, '  TO  PRINT  ALL  OF  THE  MATRICES  ENTER  1,  ELSE  ENTER  0  > 
READ(8,*)  Q 
PRINT' (///)' 

C 

C 

C  READ  IN  THE  WEIGHTING  MATRIX  DIAGONAL 
C  VALUE  FOR  EACH  MODE 

C 

C 

PRINT*, '  ENTER  THE  DIAGONAL  VALUES,  IN  MODE  INPUT  ' 

PRINT*,  *  ORDER,  FOR  THE  CONTROL  WEIGHTING  MATRIX  > ' 

READ(8,*)  (AA(I), I=1,N) 

PRINT*,'  ' 

PRINT*, (AA(I), 1=1, N) 

PRINT' (//)' 

PRINT' (///) ' 

C 

C 

C  FORMING  THE  A,B,C  AND  WEIGHTING  MATRICES 

C 

C 

CALL  FORMA(Al, D, W,NC1,NC12, IC1) 

CALL  FORMB(Bl, PHIA, NCI, NC12, NACT, IC1 ) 

CALL  FORMC(Cl, PHIS, NCI, NC12, NSEN, IC1 ) 

CALL  F0RMA(A2, D,  W, NC2, NC 22,  IC2) 

CALL  F0RMB(B2, PHIA, NC2, NC22, NACT, IC2) 

CALL  F0RMC(C2, PHIS, NC2, NC22, NSEN, IC2) 

CALL  FORMA(A3, D, W, NC3, NC32, IC3) 

CALL  F0RMB(B3,  PHIA,  NC3,  NC32, NACT, IC3) 

CALL  F0RMC(C3, PHIS, NC3, NC32, NSEN, IC3) 

CALL  F0RMA(A4, D, W, NR, NR2, IR) 

CALL  F0RMB(B4, PHIA, NR, NR2, NACT, IR) 

CALL  F0RMC(C4, PHIS, NR, NR2, NSEN, IR) 

CALL  F0RMQ(QA1,AA,NC1, IC1) 

CALL  FORMQ(QA2,AA,  NC2,  IC2) 

CALL  F0RMQ(QA3,AA,NC3,IC3) 

IF  (DEC.EQ.4)  THEN 
CALL  F0RMQ(QA4,AA,NR,IR) 

ENDIF 

C 

C 


C  PRINTING  THE  A, B,C  AND  WEIGHTING  MATRICES 


IF  (Q.EQ.l)  THEN 

PRINT*, 1  THE  CONTROLLER  1  A  MATRIX  IS  * 

CALL  PRNT(A1,NC12,NC12) 

PRINT*. '  THE  CONTROLLER  1  B  MATRIX  IS  ' 

CALL  PRNT(B1,NC12,NACT) 

PRINT*, '  THE  CONTROLLER  1  C  MATRIX  IS  ' 

CALL  PRNT(C1,NSEN,NC12) 

PRINT*, •  THE  Cl  CONTROL  WEIGHTING  MATRIX  IS  1 
CALL  PRNT (QA1,  NC12,  NC12) 

PRINT*, '  THE  CONTROLLER  2  A  MATRIX  IS  ' 

CALL  PRNT ( A2, NC22, NC22 ) 

PRINT*, '  THE  CONTROLLER  2  B  MATRIX  IS  ' 

CALL  PRNT ( B2 , NC22 , NACT ) 

PRINT*, '  THE  CONTROLLER  2  C  MATRIX  IS  ' 

CALL  PRNT (C2, NSEN, NC22 ) 

PRINT*, '  THE  C2  CONTROL  WEIGHTING  MATRIX  IS 
CALL  PRNT (QA2, NC22, NC22 ) 

PRINT*, '  THE  CONTROLLER  3  A  MATRIX  IS  ' 

CALL  PRNT (A3, NC32, NC32 ) 

PRINT*, '  THE  CONTROLLER  3  B  MATRIX  IS  ' 

CALL  PRNT ( B3 , NC32 , NACT ) 

PRINT*, '  THE  CONTROLLER  3  C  MATRIX  IS  ' 

CALL  PRNT ( C3, NSEN, NC32 ) 

PRINT*, *  THE  C3  CONTROL  WEIGHTING  MATRIX  IS 
CALL  PRNT ( QA3, NC32, NC32 ) 

IF  (NR.EQ.0)  THEN 

PRINT*, '  NO  RESIDUAL  TERMS  1 

GOTO  115 

ENDIF 

IF  (DEC.EQ.3)  THEN 

PRINT*, '  THE  A  RESIDUAL  MATRIX  IS  ' 

CALL  PRNT(A4,NR2,NR2) 

PRINT*, '  THE  B  RESIDUAL  MATRIX  IS  ' 

CALL  PRNT(B4, NR2,NACT) 

PRINT*, '  THE  C  RESIDUAL  MATRIX  IS  ' 

CALL  PRNT (C4,  NSEN,  NR2) 

ELSE 

PRINT*, '  THE  CONTROLLER  4  A  MATRIX  IS  ' 

CALL  PRNT(A4,NR2,NR2) 

PRINT*, '  THE  CONTROLLER  4  B  MATRIX  IS  ' 

CALL  PRNT(B4,NR2,NACT) 

PRINT*, '  THE  CONTROLLER  4  C  MATRIX  IS  ' 

CALL  PRNT (C4, NSEN, NR2) 

PRINT*, '  THE  C4  CONTROL  WEIGHTING  MATRIX  IS 
CALL  PRNT(QA4,NR2,NR2) 

ENDIF 


ENDIF 

115  CONTINUE 


c 

C  THIS  SECTION  GENERATES  THE  RICCATI  SOLUTIONS 
C  AND  THE  GAIN  MATRICES  OF  EACH  CONTROLLER 
C 
C 

CALL  VMULFP(B1, Bl, NC12, NACT, NC12, NCOL, NCOL, SAT, NCOL, IER) 
120  CONTINUE 
IER  =  0 
TOL  *  0.001 

PRINT*, '  THE  FOLLOWING  ARE  THE  MRIC  A+BG  1  INPUTS  ' 

PRINT' (//)' 

PRINT*, '  THE  MATRIX  A1  IS  * 

CALL  PRNT(A1,NC12,NC12) 

PRINT*, '  THE  MATRIX  SAT  (B1*B1T)  IS  ' 

CALL  PRNT(SAT,NC12,NC12) 

PRINT*, '  THE  MATRIX  QA1  IS  ' 

CALL  PRNT ( QA1 , NC12,  NC12 ) 

PRINT*, '  NC12  =  ' , NC12 
PRINT* (//) 1 

CALL  MRIC(NC12, Al, SAT, QA1, S, ABG1, TOL, IER) 

C 

C  ABG1  *  Al  +  B1G1 
C 

PRINT*, '  THE  EIGENVALUES  IF  Al  +  B1G1  ' 

PRINT' (/)' 

CALL  EIGRF(ABG1, NC12, NCOL, 0, Wl, TEN, NCOL, STOR, IER) 

DO  629  I=1,NC12 
629  PRINT*  *  ',W1(I) 

PRINT' (/)' 

PRINT*, '  THE  RICCATI  SOLUTION  OF  AC  +  BCG  tl  IS  ' 

CALL  PRNT(S, NC12, NC12) 

CALL  VMULFM(B1,  S| NC12, NACT, NC12, NCOL,  NCOL, GAIN1, NCOL, IER) 
PRINT*, '  THE  G1  GAIN  MATRIX  IS  1 
CALL  PRNT (GAIN1, NACT, NC12) 

CALL  FINDK(C1, NSEN, NC12, WORK, GAIN1, NSEN, KHAT1 ) 

PRINT*, 'THE  FEEDBACK  MATRIX  KHAT1  IS' 

CALL  PRNT(KHAT1,  NACT,  NSEN) 

CALL  ABGC(A1, Bl,  Cl,  NC12, KHAT1, NACT, NSEN, ABK1 ) 

CALL  EIGRF(ABK1, NC12, NCOL, 0, Wl, TEN, NCOL, STOR, IER) 

PRINT*, '  EIGENVALUES  OF  Al  +  B1K1C1  ARE' 

PRINT' (/)' 

DO  632  M.NC12 
632  PRINT*,'  ' , Wl ( I ) 

PRINT' (/)' 

125  CONTINUE 

CALL  VMULFP(B2, B2, NC22, NACT, NC22, NCOL, NCOL, SAT2, NCOL, IER) 
140  CONTINUE 
IER  *  0 
TOL  *  0.001 

PRINT*, '  THE  FOLLOWING  ARE  THE  MRIC  A+BG  2  INPUTS  ' 

PRINT' (//)' 

PRINT*, '  THE  MATRIX  A2  IS  ' 
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CALL  PRNT(A2,NC22,NC22) 

PRINT*,  *  THE  MATRIX  SAT2  (B2*B2T)  IS 
CALL  PRNT(SAT2, NC22,  NC22) 

PRINT*, '  THE  MATRIX  QA2  IS  ' 

CALL  PRNT ( QA2, NC22, NC22 ) 

PRINT*, '  NC22  =  ' ,NC22 
PRINT' (//)' 

CALL  MRIC(NC22, A2, SAT2, QA2, S, ABG2, TOL, IER) 

PRINT*, '  THE  EIGENVALUES  OF  A2  +  B2G2  ' 

CALL  EIGRF(ABG2,  NC22, NCOL, 0, Wl, TEN, NCOL, STOR, IER) 

PRINT' (/)' 

DO  630  I=1,NC22 

630  PRINT*,'  ' , Wl ( I ) 

PRINT' (/)' 

PRINT*, '  THE  RICCATI  SOLUTION  OF  AC  +  BCG  #2  IS  ' 

CALL  PRNT(S,NC22,NC22) 

CALL  VMULFM(B2, S, NC22, NACT, NC22, NCOL, NCOL, GAIN2, NCOL, IER) 
PRINT*, '  THE  G2  GAIN  MATRIX  IS  ' 

CALL  PRNT (GAIN2, NACT, NC22) 

CALL  FINDK(C2,NSEN,  NC22, WORK, GAIN2, NSEN, KHAT2) 

PRINT*, 'THE  FEEDBACK  MATRIX  KHAT2  IS' 

CALL  PRNT(KHAT2, NACT, NSEN) 

CALL  ABGC(A2,  B2,  C2,  NC22, KHAT2, NACT, NSEN, ABK2 ) 

CALL  EIGRF(ABK2, NC22, NCOL, 0, Wl, TEN, NCOL, STOR,  IER) 

PRINT' (/)' 

PRINT*, '  THE  EIGENVALUES  OF  A2  +  B2K2C2  ARE' 

DO  633  I=1,NC22 
633  PRINT*,'  ' , Wl ( I ) 

PRINT' (/)' 

145  CONTINUE 

CALL  VMULFP(B3, B3, NC32, NACT, NC32, NCOL, NCOL, SAT3, NCOL, IER) 
150  CONTINUE 
IER  =  0 
TOL  =  0.001 

PRINT*, '  THE  FOLLOWING  ARE  THE  MRIC  A+BG  3  INPUTS  ' 

PRINT' (//) ' 

PRINT*, '  THE  MATRIX  A3  IS  ' 

CALL  PRNT (A3, NC32, NC32 ) 

PRINT*, '  THE  MATRIX  SAT3  <B3*B3T)  IS  ' 

CALL  PRNT (SAT3,  NC32, NC32) 

PRINT*, '  THE  MATRIX  QA3  IS  ' 

CALL  PRNT (QA3, NC32, NC32 ) 

PRINT*,'  NC32  =  ',NC32 
PRINT' (//)' 

CALL  MRIC(NC32, A3, SAT3, QA3, S, ABG3, TOL, IER) 

PRINT*, '  THE  EIGENVALUES  OF  A3  +  B3G3  ' 

CALL  EIGRF(ABG3,NC32, NCOL, 0,W1, TEN, NCOL, STOR, IER) 

PRINT' (/)' 

DO  631  I=1,NC32 

631  PRINT*,'  * , Wl ( I ) 

PRINT' (/)' 

PRINT*, '  THE  RICCATI  SOLUTION  OF  AC  +  BCG  #3  IS 
CALL  PRNT(S,NC32,NC32) 
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CALL  VMULFM(B3,S, NC32,  NACT, NC32, NCOL, NCOL, GAIN3, NCOL, IER) 
PRINT*, '  THE  G3  GAIN  MATRIX  IS1 
CALL  PRNT (GAIN3, NACT, NC32) 

CALL  FINDK(C3,  NSEN,  NC32, WORK, GAIN3, NSEN, KHAT3) 

PRINT*, 'THE  FEEDBACK  MATRIX  KHAT3  IS' 

CALL  PRNT(KHAT3,  NSEN,  NSEN) 

CALL  ABGC(A3, B3, C3, NC32, KHAT3, NACT, NSEN, ABK3) 

CALL  EIGRF(ABK3, NC32, NCOL, 0, Ml, TEN, NCOL, STOR, IER) 

PRINT*, '  THE  EIGENVALUES  OF  A3  +  B3K3C3  ' 

PRINT' (/)' 

DO  634  I=1,NC32 
634  PRINT*,'  ' , W1 ( I ) 

PRINT' (/)' 

155  CONTINUE 

IF  (DEC.EQ.4)  THEN 

CALL  VMULFP(B4, B4, NR2, NACT, NR2, NCOL, NCOL, SAT4, NCOL, IER) 
160  CONTINUE 
IER  =  0 
TOL  =  0.001 

CALL  MRIC(NR2, A4, SAT4, QA4, S, ABG4, TOL, IER) 

PRINT*, '  THE  RICCATI  SOLUTION  OF  AC  +  BCG  #4  IS  ’ 

CALL  PRNT(S,NR2,NR2) 

CALL  VMULFM(B4, S, NR2, NACT, NR2, NCOL, NCOL, GAIN4, NCOL, IER) 
PRINT*, '  THE  G4  GAIN  MATRIX  IS  ' 

CALL  PRNT ( GAIN4, NACT , NR2 ) 

CALL  FINDK(C4,  NSEN,  NR2, WORK, GAIN4, NSEN, KHAT4) 

PRINT*, 'THE  FEEDBACK  MATRIX  KHAT4  IS' 

CALL  PRNT(KHAT4, NSEN, NSEN) 

ENDIF 

165  CONTINUE 

CALCULATE  SUMMATION  OF  KHATS 

CALL  MADD( KHAT1 ,  KHAT2 ,  NACT ,  NACT ,  KCC ) 

CALL  MADD(KCC,KHAT3,  NACT, NACT,  BCG) 

IF(DEC.EQ.3)THEN 
DO  166  1=1, NACT 
DO  166  0=1, NACT 

166  KCC(I, J)=BCG(I, J) 

ELSE 

CALL  MADD(BCG,KHAT4, NACT, NACT, KCC) 

ENDIF 

SUMMATION  OF  KHAT  IS  NOW  IN  KCC 

THIS  SECTION  GENERATES  THE  BLOCK  SEGMENTS 
OF  MAOM  ANO  PUTS  THEM  INTO  THE  MAOM  MATRIX 

THE  FOUR  CONTROLLER  MATRIX  WILL  CONTAIN 
RESIDUAL  TERMS  (SEE  DIAGRAM  BELOW). 

THE  FOUR  CONTROLLER  MATRIX  DOES  NOT  IN¬ 
CLUDE  RESIDUALS  (YET). 


c  THE  FOUR  CONTROLLER  MATRIX  (MAUM)  WITH 
C  RESIDUAL  TERMS  WILL  LOOK  LIKE: 

C 

c  **** 


c 

* 

* 

c 

* 

A1+B1KIC1 

B1KIC2 

B1KIC3 

B1KIC4 

B1KICR 

c 

* 

c 

B2KIC1 

A2+B2KIC2 

B2KIC3 

B2KIC4 

B2KICR 

* 

c 

» 

* 

c 

* 

B3KIC1 

B3KIC2 

A3+B3KIC3 

B3KIC4 

B3KICR 

* 

c 

* 

* 

c 

* 

B4KIC1 

84KIC2 

B4KIC3 

A4+B4KIC4 

B4KICR 

* 

c 

* 

* 

c 

* 

BRKIC1 

BRKIC2 

BRKIC3 

BRKIC4 

AR+BRKICR 

* 

c 

* 

* 

c 


c 

K=NC12 
L=NC22  +  K 
P5=NC32  +  L 

MM  =  NC12  +  NC22  +  NC32  +  NR2 
C 

DO  200  1=1, MM 
DO  200  J=1,MM 

200  MAOM(I.O)  =  0.0 

CALL  ABGC( Al, Bl, Cl , NC12 , KCC, NACT, NSEN, ABG1 ) 
DO  201  1=1, NC12 
DO  201  J-1.MC12 

201  MAOM(I,J)  =  ABQ1(I,J) 

CALL  ABGC(A2  82, C2, NC22, KCC, NACT, NSEN, ABG2) 
DO  202  I=1,NC22 
DO  202  J=1,NC22 

202  MAJM(I+K, J+K)  =  ABG2(I,0) 

CALL  ABGC(A3, B3, C3, NC32, KCC, NACT, NSEN, ABG3) 
DO  203  I=1,NC32 
DO  203  J=1,NC32 

203  MA0M(I+L,0+L)  =  ABG3(I,J) 

CALL  MMUL(B1, KCC, NC12, NACT, NACT, TEN) 

CALL  MMUL(TEN,C2,NC12, NACT, NC22, BCG) 

DO  208  I=1,NC12 
DO  208  0=1,NC22 

208  MAJM(I,J*K)  =  BCG(I, J) 

CALL  MMUL(B1, KCC,  NC12,  NACT, NACT, TEN) 

CALL  MMUL(TEN,C3,NC12, NACT, NC32, BCG) 

DO  209  1=1, NC12 
DO  209  J=1,NC32 

209  MA0M(I,0+L)  =  BCG(I, J) 

CALL  MMUL(B2, KCC, NC22, NACT, NACT, TEN) 

CALL  MMUL(TEN, Cl, NC22, NACT, NC12, BCG) 

00  210  I=1,NC22 
DO  210  0=1,NC12 

210  MA0M(I+K,3)  =  BCG(I,0) 
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CALL  MMUL(B2, KCC, NC22, NACT, NACT, TEN) 

CALL  MMUL(TEN,  C3,  NC22, NACT, NC32, BCG) 

00  212  I=1,NC22 
00  212  J=1 , NC32 

212  MA0M(I+K, 0+L)  =  BCG(I,3) 

CALL  MMUL(B3, KCC,  NC32, NACT, NACT, TEN) 

CALL  MMUL(TEN,C1,NC32,NACT,NC12,BCG) 

DO  213  I=1,NC32 
DO  213  0=1,NC12 

213  MA0M(I+L,0)  =  BCG(I,0) 

CALL  MMUL(B3, KCC, NC32, NACT, NACT, TEN) 

CALL  MMUL(TEN,C2,NC32,  NACT, NC22, BCG) 

DO  214  I=1,NC32 
DO  214  0=1,NC22 

214  MA0M(I+L,0+K)  =  BCG(I,J) 

CALL  MMUL(B4, KCC, NR2, NACT, NACT, TEN) 

CALL  MMUL(TEN,C1,NR2,NACT,NC12,BC6) 

DO  222  1*1, NR2 
DO  222  J=1,NC12 

222  MA0M(I+P5,  J)  =*  BCG(I,J) 

CALL  MMUL(B4, KCC, NR2, NACT, NACT, TEN) 

CALL  MMUL(TEN,C2,NR2, NACT, NC22, BCG) 

DO  223  1=1, NR2 
DO  223  J=l, NC22 

223  MAJM(I+P5, J+K)  =  BCG(I.J) 

CALL  MMUL(B4,  KCC, NR2, NACT, NACT, TEN) 

CALL  MMUL(TEN, C3,  NR2, NACT, NC32,  BCG) 

DO  224  1=1, NR2 
DO  224  J=1,NC32 

224  MAJM(I+P5,J+L)  =  BCG(I,0) 

C 

C 

IF  (DEC.EQ.4)  THEN 
C 
C 

CALL  ABGC(A4, B4, C4, NR2, KCC, NACT, NSEN, ABG4) 
DO  230  1=1, NR2 
DO  230  0=1, NR2 

230  MA0M(I+P5,O*P5)  =  ABG4(I,0) 

CALL  MMUL(B1, KCC, NC12, NACT, NACT, TEN) 

CALL  MMUL(TEN,C4,NC12, NACT, NR2, BCG) 

DO  232  I=1,NC12 
DO  232  0=1, NR2 

232  MAOM(I,0+P5)  =  BCG(I,0) 

CALL  MMUL(B2, KCC,  NC22, NACT, NACT, TEN) 

CALL  MMUL(TEN,C4,NC22, NACT, NR2, BCG) 

DO  233  I=1,NC22 
DO  233  0=1, NR2 

233  MA0M(I+K,0+P5)  =  BCG(I,0) 

CALL  l«UL(B3, KCC, NC32, NACT, NACT, TEN) 

CALL  MMUL(TEN,C4,NC32, NACT, NR2, BCG) 

DO  234  I=1,NC32 
DO  234  0=1, NR2 
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234  MAJM(I+L, J+P5)  =  BCG(I, J) 

C 

ELSE 

C 

CALL  F0RMA(A4,D, W, NR,  NR2, IR) 

CALL  ABGC(A4, B4, C4, NR2, KCC, NACT, NSEN, ABG4) 
DO  250  1=1, NR2 
DO  250  J=1,NR2 

250  MAJM(I+P5, J+P5)  =  ABQ4(I, J) 

ENDIF 


NOW  WE  HAVE  THE  MAJM  MATRIX 


IF  (DEC.EQ.4)  THEN 

PRINT*, 1  THE  FOUR  CONTROLLER  MAJOR  MATRIX  ' 

ELSE 

PRINT*, '  THE  THREE  CONTROLLER  MAJOR  MATRIX  ' 
ENDIF 

CALL  PRNTXL(MAJM,MM,MM) 


NEXT,  FORM  E  TO  THE  AT 


PRINT*, 1  FOR  THE  TIME  RESPONSE  AND  THE  EIGENVALUE' 

*  , 1  ANALYSIS  ENTER  1  1 

PRINT*, '  FOR  ONLY  THE  EIGENVALUE  ANALY5IS  ENTER  2  >  ' 

READ(8,*)  SKIP 
IF  (SKIP.EQ.2)  THEN 
PRINT*,'  ' 

PRINT*, '  YOU  CHOSE  "2”  SO  WE  WILL  FORGET  THE  TIME  RESPONSE  ' 

GOTO  300 

ENDIF 

PRINT*, '  YOU  WANT  BOTH,  EH?  ' 

PRINT*, '  ' 

PRINT*, '  0. K.  HERE  THEY  ARE!  ' 

PRINT*,'  ' 

PRINT*, '  THE  TIME  RESPONSE  IS  FIRST. . .  ’ 

DT  =  0. 1 
TOL  =  0.001 

CALL  MEXP(MM, MAJM, DT, EAT2) 

PRINT*,'  MEXPOUT  ' 

PRINT' (//)' 

PRINT*, '  THE  SOLUTION  EAT2  IS  ' 

CALL  PRNTXL ( EAT2, MM, MM) 


EAT2  IS  NOW  THE  SOLUTION  TO  E  TO  THE  AT 


CALL  F0RMX2(X0, INIT) 
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m  «a™E(»«,  MM.  or,  XI,  XO,  MODE,  EAT,  WORK,  DEC ) 


C  EIGENVALUE  ANALYSIS  SECTION 

C 

C 

PRINT '(///)* 

PRINT*, 1  OVERALL  SYSTEM  EIGENVALUES  • 

CALLtEIGRP (MAJM, MM,  NDA, 0, Z, TEN, NCOL,  WORK,  IER) 

DO  400  1=1, MM 

400  PRINT*,'  'Z(I) 

PRINT' (//)' 

C 

PRINT*, '  EIGENVALUES  OF  A1  +  B1KIC1 ' 

PRINT' (/V(ABG1,  NC12’ NC0L’ 0’ W1,  TEN’ NC0L’  ST0R’ IER) 

DO  401  1=1, NC12 

401  PRINT*,'  \W1(I) 

PRINT' (//)' 

C 

PRINT*, '  EIGENVALUES  OF  A2  +  B2KIC2' 

PRINT  (/ V ^(ABG2’ NC22’ NC0L’ W1, TEN'  NC0L’  ST0R» IER) 

DO  403  I=1,NC22 
403  PRINT*,'  ' , W1 ( I ) 

PRINT* (//)' 

PRINT*, '  EIGENVALUES  OF  A3  +  B3KIC3' 

PRINT^( / V <AB63’ NC32' NC0L’ 0>  W1>  TEN’ NC°L' 5T0R? IER) 

DO  405  I=1,NC32 
405  PRINT*,'  '  W1 ( I ) 

PRINT' (//)' 

IF  (NR.EQ.0)  THEN 

PRINT*, '  NO  RESIDUAL  TERM  EIGENVALUES  ' 

PRINT' (/)' 

GOTO  410 
ENDIF 

IF  (DEC. EQ. 4)  THEN 
C 

PRINT*, '  EIGENVALUES  OF  A4  +  B4KIC4' 

CALL  EIGRF (ABG4, NR2, NCOL, 0, Wl, TEN, NCOL,  STOR,  IER) 

PRINT  ( / ) 

DO  407  1=1, NR2 

407  PRINT*,'  ' , Wl ( I ) 

PRINT' (//) ' 

ELSE 

C 

.  EIGENVALUES  OF  THE  A  RESIDUAL  MATRIX 

CALL  EIGRF (A4, NR2, NCOL, 0, Wl, TEN, NCOL, STOR, IER) 

PRINT ' ( / ) 1 
DO  409  1=1, NR2 
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409  PRINT*, 

ENDIF 
PRINT '(///)' 

C 

410  CONTINUE 


THIS  SECTION  FORMS  THE  TRANSFORMATION  MATRICES. 


AFTER  THE  TRANSFORMATION  IS  COMPLETE, 

THE  THREE  OR  FOUR  CONTROLLER  MAJM  (WITHOUT  RESIDUALS) 
WILL  LOOK  LIKE: 


* 

*  (A1+B1K1C1) 

0 

0 

* 

0 

*  (B2K1C1) 

(A2+B2K2C2) 

0 

0  * 
* 

w 

*  (B3K1C1+B3K2C1) 

(B3K2C2) 

(A3+B3K3C3) 

0  * 

*  ( B4K1C1+B4K2C1+B4K3C1 ) 

(B4K2C2+B4K3C2) 

( B4K3C3 ) 

(A4+B4K4C4)  * 
* 

WHERE  THE  NON-ZERO  TERMS  INCLUDE  THE 
TRANSFORMATION  MATRICES. 

(THE  RESIDUALS  ARE  NOT  ZEROED  SINCE  THEY 
ARE  NOT  ADDRESSED  IN  THIS  PROGRAM) 


GENERATE  THE  TRANSFORMATION  MATRICES 

GAMMA1,  GAMMA2,  GAMMA3  (GAMMA4  =  IDENTITY) 

CALL  TFR(CT,C2,NSENfNC22,l,2) 

DO  500  1=1, NC2 
DO  500  0=1,NSEN 

500  V(I,0)  =  CT(I+NC2, 0) 

CALL  TFR ( CT, C3, NSEN, NC32, 1,2) 

DO  501  1=1, NC3 
DO  501  J=1 , NSEN 

501  V(I+NC2,0)  =  CT(I+NC3, J) 

IF  (DEC.EQ.4)  THEN 

CALL  TFR(CT,C4, NSEN, NR2, 1,2) 

DO  502  1=1, NR 
DO  502  0=1, NSEN 

502  V(I+NC2+NC3,0)  =  CT(I+NR,J) 

NRV  =  NC2  +  NC3  +  NR 
PRINT*,'  V  (C2/C3/C4)  IS  ' 
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ELSE 

PRINT*,’  V  (C2/C3)  IS  1 
NRV  =  NC2  +  NC3 
ENDIF 

CALL  PRNT(V, NRV, NSEN) 

CALL  LSVDF( V, NCOL, NRV, NSEN, TEN, NCOL, -1, SING, STOR, IER ) 
PRINT*,'  ' 

PRINT*, '  V  OUT  OF  LSVOF  IS  ’ 

CALL  PRNT(V, NSEN, NSEN) 

PI  =  NSEN  -  NRV 
IF  (Pl.LT.l)  THEN 
DO  503  1=1, NSEN 

503  GAMMA1(I, 1)  =  V(I,NSEN) 

PI  =  1 

ELSE 

DO  504  1=1, NSEN 
DO  504  0=1, PI 

504  GAMMA1(1,J)  =  V(I,0+NRV) 

ENDIF 

PRINT*, '  TRANSFORMATION  MATRIX  GAMMA1  ' 

CALL  PRNT ( GAMMA 1, NSEN, PI ) 

CHECK  TO  SEE  THAT  GAMMA1  IS  ORTHOGONAL  TO  BOTH  C2  AND  C3 

NOTE:  AKC  IN  THIS  SECTION  IS  OUST  A  WORK  AREA  TO  TEST 

THE  ORTHOGONALITY  OF  CT  *  TR.  IN  ALL  CASES  IT 
SHOULD  BE  A  BLOCK  ZERO  MATRIX. 

CALL  TFR(CT, C2, NSEN, NC22, 1,2) 

CALL  MMUL ( CT, GAMMA1 , NC22, NSEN, PI , AKC ) 

PRINT*, '  C2T  *  GAMMA1  1 
CALL  PRNT(AKC,NC22,P1) 

CALL  TFR(CT, C3, NSEN, NC32, 1,2) 

CALL  MMUL ( CT, GAMMA 1 , NC32, NSEN,  PI , AKC ) 

PRINT*, '  C3T  *  GAMMA1  1 
CALL  PRNT(AKC,NC32,P1) 

IF  (DEC.EQ.4)  THEN 

CALL  TFR(CT,C4, NSEN, NR2, 1,2) 

CALL  MMUL ( CT, GAMMA1 , NR2 , NSEN, PI , AKC ) 

PRINT*, 1  C4T  *  GAMMA1  ' 

CALL  PRNT (AKC, NR2, PI) 

ENDIF 

IF  (DEC.EQ.4)  THEN 

PRINT*, '  C234  SINGULAR  VALUES  ' 

- ELSE . 

PRINT*, ’  C23  SINGULAR  VALUES  ' 

ENDIF 

CALL  PRNT (SING, NRV, 1) 

CALL  TFR(TRT,GAMMA1,NSEN,P1, 1,2) 

CALL  MMUL(TRT,GAMMA1,P1,NSEN,P1,RK) 

CALL  GMINV(P1,P1,RK,RK1, J,TAPE) 

CALL  TFR(CT,C1,NSEN, NC12, 1,2) 
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CALL  MMUL ( TRT, Cl, PI, NSEN, NC12, AKC ) 

CALL  MMUL ( CT, GAMMA 1 ,  NC12, NSEN, PI , K0B1 ) 
CALL  MMUL(K0B1,RK1,NC12,P1,P1,ST0R) 
CALL  MMUL(STOR, AKC, NC12, PI, NC12, CTCC1 ) 


CALL  TFR(CT,C3, NSEN, NC32, 1,2) 

DO  505  1=1, NC3 
DO  505  J=1,NSEN 

505  V(I,0)  =  CT(I+NC3,0) 

IF  (DEC.EQ.4)  THEN 

CALL  TFR(CT, C4, NSEN, NR2, 1,2) 

DO  506  1=1, NR 
DO  506  J=l, NSEN 

506  V(I+NC3, J)  =  CT(I+NR, 0) 

NRV  =  NC3  +  NR 
PRINT*, '  V  (C3/C4)  IS  ' 

ELSE 

NRV  =  NC3 

PRINT*,  1  V  (C3)  IS  * 

ENDIF 

CALL  PRNT(V, NRV, NSEN) 

CALL  LSVDF(V, NCOL, NRV, NSEN,  TEN, NCOL,  -1, SING,  STOR,  IER) 
PRINT*,'  ' 

PRINT*, '  V  OUT  OF  LSVDF  IS  ' 

CALL  PRNT(V, NSEN, NSEN) 

P2  =  NSEN  -  NRV 
IF  (P2.LT.1)  THEN 
DO  507  1=1, NSEN 

507  GAMMA2 (1,1)  =  V(I,NSEN) 

P2  =  1 

ELSE 

DO  508  1=1, NSEN 
DO  508  0=1, P2 

508  GAMMA2(I, 0)  =  V(I,0*NRV) 

ENDIF 

PRINT*, 1  TRANSFORMATION  MATRIX  GAMMA2  ' 

CALL  PRNT(GAMMA2,NSEN,P2) 

C 

C  CHECK  TO  SEE  THAT  GAMMA2  IS  ORTHOGONAL  TO  C3 
C 

CALL  TFR(CT,C3, NSEN, NC32, 1,2) 

CALL  MMUL(CT, GAMMA2.NC32, NSEN,  P2, AKC) 

PRINT*, '  C3T  *  GAMMA2  ' 

CALL  PRNT(AKC,NC32, P2) 

IF  (DEC.EQ.4)  THEN 

CALL  TFR (CT,C4, NSEN, NR2, 1,2) 

CALL  MMUL(CT,GAMMA2,NR2, NSEN, P2, AKC) 

PRINT*, '  C4T  *  GAMMA2  ' 

CALL  PRNT(AKC,NR2,P2) 


86 


>  r~>  o 


ENDIF 

IF  (DEC. EQ.4)  THEN 

PRINT*, •  C34  SINGULAR  VALUES  ' 

ELSE 

PRINT*, '  C3  SINGULAR  VALUES  ' 

ENDIF 

CALL  PRNT(SING,NRV,1) 

CALL  TFR(TRT, GAMMA2, NSEN, P2, 1, 2) 

CALL  MMUL(TRT,GAMMA2,P2,NSEN,  P2,RK) 
CALL  GMINV(P2,P2, RK, RK2, 0, TAPE) 

CALL  TFR (CT,C2, NSEN, NC22, 1,2) 

CALL  MMUL(TRT, C2, P2, NSEN, NC22, AKC) 

CALL  MMUL ( CT , GAMMA2, NC22, NSEN, P2, K0B2 ) 
CALL  MMUL ( KOB2 , RK2 ,NC22,P2,P2,STOR) 
CALL  MMUL (STOR, AKC, NC22, P2.NC22, CTCC2) 


IF  (DEC. EQ.4)  THEN 

CALL  TFR(CT,C4, NSEN, NR2, 1,2) 

DO  509  1*1, NR 
DO  509  J=1,NSEN 

509  V(I,0)  =  CT{I+NR, J) 

PRINT*,'  V  (C4)  IS  ' 

CALL  PRNT(V,NR,NSEN) 

CALL  LSVDF( V, NCOL, NR, NSEN, TEN, NCOL, -1, SING, STOR, IER) 
PRINT*, '  V  OUT  OF  LSVDF  IS  1 
CALL  PRNT(V, NSEN, NSEN) 

P3  =  NSEN  -  NR 
IF  (P3.LT.1)  THEN 
DO  510  1=1, NSEN 

510  GAMMA3(I, 1)  =  V(I,NSEN) 

P3  =  1 

ELSE 

DO  511  1=1, NSEN 
DO  511  J=1,P3 

511  GAMMA3(I,0)  =  V(I,J+NR) 

ENDIF 

PRINT*, '  TRANSFORMATION  MATRIX  GAMMA3  ' 

CALL  PRNT ( GAMMA3, NSEN,  P3 ) 

CHECK  TO  SEE  THAT  GAMMA3  IS  ORTHOGONAL  TO  C4 

CALL  MMUL(CT, GAMMA3, NR2, NSEN, P3, AKC) 

PRINT*, '  C4T  *  GAMMA3  ' 

CALL  PRNT(AKC,NR2,P3) 

PRINT*, '  C4  SINGULAR  VALUES  ' 
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CALL  PRNT(SING,  NR, 1) 

CALL  TFR(TRT,GAMMA3, NSEN, P3, 1,2) 

CALL  MMUL ( TRT, GAMMA3, P3, NSEN,  P3,  RK ) 
CALL  GMINV ( P3, P3, RK, RK3,  J,  TAPE) 

CALL  TFR(CT,C3, NSEN, NC32, 1,2) 

CALL  MMUL(TRT,C3,P3,NSEN,NC32,AKC) 

CALL  MMUL ( CT, GAMMA3 , NC32 , NSEN, P3 , K0B3 ) 
CALL  MMUL  (K0B3,  RIG,  NC32, P3,P3, STOR) 
CALL  MMUL (STOR, AKC, NC32, P3, NC32, CTCC3 ) 
ENOIF 

DO  411  1=1, NSEN 
DO  411  J=1,P1 

411  GTl(sJ,  I)=GAMMA1(I,  J) 

DO  412  1=1, NSEN 

DO  412  3=1, P2 

412  GT2(J, I)=GAMMA2(I, J) 

PRINT*, '  GAMMA1  TRANSPOSE  1 
CALL  PRNT(GT1,P1,NSEN) 

PRINT*, '  GAMMA2  TRANSPOSE  ' 

CALL  PRNT(GT2,P2,NSEN) 

IF(DEC. EQ.4)  THEN 

DO  413  1=1, NSEN 
DO  413  J=1,P3 

413  GT3(0, I)=GAMMA3(I, J) 

PRINT*, '  GAMMA3  TRANSPOSE  ' 

CALL  PRNT(GT3,P3,NSEN) 

ENDIF 

CSTAR  =  GAMMA  *  C 

CALL  MMUL(GT1, Cl, PI, NSEN, NC12, CSTAR1 ) 
PRINT*,'  ( GT1 ) (Cl )  ' 

CALL  PRNT (CSTAR1, PI, NC12) 

CALL  MMUL(GT2, C2, P2, NSEN, NC22, CSTAR2) 
PRINT*,'  (GT2)(C2)  ' 

CALL  PRNT(CSTAR2, P2,  NC22 ) 

IF(DEC. EQ.4)  THEN 

CALL  MMUL ( GT3, C3, P3, NSEN, NC32, CSTAR3 ) 
PRINT*,'  ( GT  3 ) ( C3 )  ' 

CALL  PRNT ( CSTAR3, P3, NC32 ) 

ENDIF 


TRANSFORMATION  MATRICES  T2,  T3,  T4  (T1  =  IDENTITY) 


DO  512  1=1, NCI 
DO  512  J=1,NACT 
512  V(I,J)  =  B1(I+NC1, J) 

PRINT*, '  V  (Bl)  IS  ' 

CALL  PRNT (V, NCI, NACT) 

CALL  LSVDF ( V , NCOL , NC 1 , NACT,  TEN, NCOL , - 1 , SING, STOR , I ER ) 
PRINT*, '  V  OUT  OF  LSVDF  IS  ' 
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CALL  PRNT ( V, NACT, NACT) 

E2  =  NACT  -  NCI 
IF  (E2.LT.1)  THEN 
00  513  1=1, NACT 

513  T2(I, 1)  =  V ( I , NACT ) 

E2  =  1 

ELSE 

00  514  1=1, NACT 
DO  514  0=1, E2 

514  T2(I,0)  =  V(I, 0+NC1) 

ENDIF 

PRINT*, '  TRANSFORMATION  MATRIX  T2  1 
CALL  PRNT(T2,NACT,E2) 

CHECK  TO  SEE  THAT  T2  IS  ORTHOGONAL  TO  B1 

NOTE:  IN  THIS  SECTION,  BCG  IS  THE  WORK  AREA 

FOR  B  *  T.  IN  ALL  CASES  THESE  SHOULD 
BE  BLOCK  ZERO  MATRICES. 


CALL  MMUL(B1,T2,NC12, NACT,  E2,  BCG) 

PRINT*, '  B1  *  T2  ’ 

CALL  PRNT(BCG, NC12.E2) 

PRINT*, '  B1  SINGULAR  VALUES  1 
CALL  PRNT (SING, NCI, 1) 

CALL  VMULFM(T2, T2, NACT, E2, E2, NACT, NACT , RK, NACT, IER ) 

CALL  GMINV(E2,E2,RK, RG2, J, TAPE) 

CALL  MMUL(B2,T2,NC22,NACT,E2,K0B2) 

CALL  MMUL(K0B2, RG2, NC22, E2, E2, SAT2) 

CALL  VMULFP(SAT2, T2, NC22,  E2, NACT, NCOL, NACT, K0B2, NCOL, IER) 
CALL  VMULFP ( K0B2 , B2 , NC22 , NACT ,  NC22 , NCOL , NCOL , SAT 2 ,  NCOL , I ER ) 


THIS  SAT2  WILL  BE  SUBSTITUTED  BACK  INTO  MR1C 
SYSTEM  2  FOR  A  NEW  G2. 


TRANSFORMATION  MATRIX  T3 


DO  515  1=1, NCI 
DO  515  0=1, NACT 

515  V(I,0)  =  B1(I+NC1,0) 

DO  516  1=1, NC2 

DO  516  0=1, NACT 

516  V(I+NC1, 0)  =  B2(I+NC2, 0) 

V  IS  NOW  AN  NC1+NC2  BY  NACT  MATRIX 
PRINT*, '  V  (B1/B2)  IS  ' 
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NRV  =  NCI  +  NC2 
CALL  PRNT(V, NRV,  NACT) 

CALL  LSVDF( V, NCOL, NRV, NACT,  TEN, NCOL, -1 , SING,  STOP,  IER ) 

PRINT*, '  V  OUT  OF  L5VDF  IS  1 
CALL  PRNT(V, NACT, NACT) 

E3  =  NACT  -  NRV 
IF  (E3.LT.1)  THEN 
DO  517  I  =1,NACT 

517  T3( 1,1)  =  V(I,NACT) 

E3  =  1 

ELSE 

DO  518  1=1, NACT 
DO  518  0=1, E3 

518  T3(I,J)  =  V ( I , O+NRV ) 

END  IF 

PRINT*, '  TRANSFORMATION  MATRIX  T3  1 
CALL  PRNT(T3,NACT,E3) 

C 

C  CHECK  TO  SEE  THAT  T3  IS  ORTHOGONAL  TO  B1  AND  B2 
C 

CALL  MMUL(B1, 13, NC12,NACT, E3, BCG) 

PRINT*, '  B1  *  T3  ' 

CALL  PRNT (BCG, NC12, E3) 

CALL  MMUL(B2,T3,NC22,  NACT,  E3,  BCG) 

PRINT*, '  B2  *  T3  ' 

CALL  PRNT(BCG, NC22.E3) 

C 

PRINT*, '  B12  SINGULAR  VALUES  ' 

CALL  PRNT (SING, NRV, 1) 

C 

CALL  VMULFM(T3, T3, NACT, E3,  E3, NACT, NACT, RK, NACT,  IER) 

CALL  GMINV ( E3, E3, RK, RG3 , 0, TAPE ) 

CALL  MMUL(B3,T3,NC32,NACT,E3,K0B3) 

CALL  MMUL(K0B3, RG3, NC32, E3,  E3,  SAT3) 

CALL  VMULFP(SAT3, T3, NC32, E3, NACT, NCOL, NACT, K0B3, NCOL, IER) 
CALL  VMULFP(K0B3, B3, NC32, NACT, NC32, NCOL, NCOL, SAT3, NCOL, IER) 
C 
C 

IF  (DEC.EQ.4)  THEN 
DO  519  1=1, NCI 
DO  519  0=1, NACT 

519  V(I,0)  =  B1(I+NC1,0) 

DO  520  1=1, NC2 

DO  520  0=1, NACT 

520  V(I+NC1,0)  =  B2(I+NC2,0) 

DO  521  1=1, NC3 

DO  521  0=1, NACT 

521  V(I+NC1+NC2,0)  =  B3(I+NC3,0) 

V  IS  NOW  A  (NC1+NC2+NC3)  BY  NACT  MATRIX 

PRINT*,'  V(B1/B2/B3)  IS  ' 

NRV  =  NCI  +  NC2  +  NC3 
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CALL  PRNT (V, NRV, NACT) 

CALL  LSVDF( V, NCOL, NRV, NACT,  TEN,  NCOL, -1 . SING,  STOR,  IER ) 
PRINT*, 1  V  OUT  OF  LSVDF  IS  ' 

CALL  PRNT (V, NACT, NACT) 

E4  =  NACT  -  NRV 
IF  (E4.LT.1)  THEN 
DO  522  1=1,  NACT 

522  T4(l, 1)  =  V(I,NACT) 

E4  =  1 

ELSE 

DO  523  1=1, NACT 
DO  523  J=1,E4 

523  T4(I, J)  =  V(I, J+NRV) 

ENDIF 

PRINT*, '  TRANSFORMATION  MATRIX  T4  ' 

CALL  PRNT(T4, NACT, E4) 

C 

C  CHECK  TO  SEE  THAT  T4  IS  ORTHOGONAL  TO  B1,B2,B3 

C 

CALL  MMUL(B1,T4,NC12, NACT, E4, BCG) 

PRINT*, •  B1  *  T4  ' 

CALL  PRNT(BCG, NC12.E4) 

CALL  MMUL(B2,T4,NC22, NACT,  E4, BCG) 

PRINT*, '  B2  *  T4  1 
CALL  PRNT (BCG, NC22, E4) 

CALL  MMUL(B3,T4,NC32, NACT,  E4,  BCG) 

PRINT*, '  B3  *  T4  ' 

CALL  PRNT(BCG,NC32,E4) 

C 

PRINT*, '  B123  SINGULAR  VALUES  ' 

CALL  PRNT (SING, NRV, 1) 

C 

CALL  VMULFM(T4, T4, NACT, E4, E4, NACT, NACT, RK, NACT, IER) 

CALL  GMINV(E4, E4, RK, RG4, 0, TAPE) 

CALL  MMUL(B4, T4, NR2, NACT, E4,  K0B4) 

CALL  MMUL(K0B4, RG4, NR2, E4,  E4,  SAT4) 

CALL  VMULFP(SAT4, T4, NR2, E4, NACT, NCOL, NACT, K0B4, NCOL,  IER ) 
CALL  VMULFP(K0B4, B4, NR2, NACT, NR2, NCOL, NCOL, SAT4, NCOL, IER) 
ENDIF 


GSTAR1  =  GAIN1 

CALCULATE  UT,  V,  QPLUS  WHERE  CSTAR  =  U*Q*VT  AND 

CSTAR+  =  V*QPLUS*UT 
KSTAR  =  GSTAR*CSTAR+ 

CALL  FINDK(CSTAR1, PI, NC12, WORK, GAIN1 , NACT, KSTAR 1 ) 
PRINT*, '  KSTAR 1  ' 

CALL  PRNT ( KSTAR1 , NACT , PI ) 
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CALL  MMUL{B2, T2, NC22,  NACT,  £2,  BSTAR2) 

ast  SS® ^  ™ 

PRINT*'(  /)  *THE  £IQENVALUES  0F  A2  +  BS1AR2QSTAR2  ' 

CALL  EIGRF(ABG2, NC22, NCOL,  0, Wl, TEN, NCOL, STOR,  IER) 

DO  714  I=1,NC22 
PRINT*,'  ' , Wl ( I ) 

PRINT'!/)' 

CALL  VMULFM(BSTAR2, S, NC22, E2,  NC22, NCOL, NCOL, GSTAR2, NCOL, IER) 
PRINT*, '  THE  GSTAR2  MATRIX  1 


CALL  PRNT ( GSTAR2, E2, NC22 ) 

PRINT' (/)' 

CALL  FINDK{CSTAR2, P2, NC22, WORK, GSTAR2, E2, KSTAR2) 
PRINT*, '  KSTAR2  * 

CALL  PRNT (KSTAR2, E2, P2 ) 


CALL  MMUL(B3, T3, NC32, NACT,  E3,  BSTAR3) 

CALL  VMULFP(BSTAR3, BSTAR3, NC32,  E3, NC32, NCOL, NCOL, SAT3, NCOL, IER) 
CALL  MRIC(NC32,A3,SAT3,QA3,S,ABG3,TOL, IER) 

PRINT*, '  THE  EIGENVALUES  OF  A3  +  BSTAR3GSTAR3  ’ 

PRINT'(/)' 

CALL  EIGRF(ABG3, NC32,  NCOL, 0,  Wl, TEN, NCOL, STOR, IER) 

DO  722  1=1, NC32 
PRINT*,'  ',W1(I) 

PRINT'!/)' 

CALL  VMULFM(BSTAR3, S, NC32, E3, NC32, NCOL, NCOL, GSTAR3, NCOL, IER) 
PRINT*, '  THE  GSTAR3  MATRIX  ' 

CALL  PRNT ( GST AR3, E3, NC32) 

PRINT'!/)' 

IF(DEC. EQ. 3)  THEN 


CALL  FINDK(C3, NSEN, NC32, WORK, GSTAR3, E3, KSTAR3) 

ELSE 

CALL  FINDK(CSTAR3, P3, NC32, WORK, GSTAR3, E3.KSTAR3) 

CALL  MMUL(B4, T4, NR2, NACT, E4, BSTAR4) 

CALL  VMULFP(BSTAR4,BSTAR4,NR2,E4,NR2, NCOL, NCOL, 5AT4, NCOL, IER) 
CALL  MRIC(NR2, A4, SAT4, QA4, S, ABG4, TOL, IER) 

CALL  VMULFM(BSTAR4, S, NR2, E4, NR2, NCOL, NCOL, GSTAR4, NCOL, IER) 
PRINT*, '  EIGENVALUES  OF  A4  +  BSTAR4GSTAR4  ' 

PRINT'!/)' 

CALL  EIGRF(ABG4, NR2, NCOL, 0, Wl, TEN, NCOL, STOR,  IER) 

DO  420  1=1, NR2 
PRINT*,'  \W1(I) 

PRINT'!/}' 

CALL  FINDK(C4, NSEN, NR2,  WORK, GSTAR4, E4, KSTAR4) 

ENOIF 


FORM  TKG  =  T  *  KSTAR  *  GAMMA 

CALL  MMUL ( KSTAR1, GT1 , NACT, PI , NSEN, TKG1 ) 
PRINT*,'  (KSTAR1)(GT1)  ' 

CALL  PRN1(TKG1, NACT, NSEN) 


CALL  MMUL i T2, KSTAR2, NACT,  £2, P2, BCG) 

CALL  MMUL (BCG, GT2, NACT , P2, NSEN, TKG2) 

PRINT*,1  (T2)(KSTAR2)(GT2)  ' 

CALL  PRNT(TKG2, NACT, NSEN) 

IF (DEC. EQ. 3)  THEN 

CALL  MMUL(T3, KSTAR3,  NACT,  E3,  NSEN, TKG3) 

PRINT*,'  (73)(KSTAR3)  1 

CALL  PRNT(TKG3, NACT, NSEN) 

ELSE 

CALL  MMUL(T3,KSTAR3, NACT, E3,P3, BCG) 

CALL  MMUL (BCG, GT3, NACT, P3, NSEN, TKG3) 

CALL  MMUL(T4, KSTAR4, NACT, E4, NSEN, TKG4) 

PRINT*,'  (T3)(KSTAR3)(GT3)  ' 

CALL  PRNT(TKG3, NACT, NSEN) 

PRINT*,'  (T4)(KSTAR4)  ' 

CALL  PRNT(TKG4, NACT, NSEN) 

ENDIF 

C 

C  FORM  (B*)(K*)(C*)  =  (B) (TKG) (C)  AND  A  +  (B*)(K*)(C*) 
C 

CALL  MMUL(B1, KSTAR1, NC12, NACT, PI, BCG) 

CALL  MMUL (BCG, CSTAR1, NC12, PI, NC12, KCC) 

PRINT*,'  ( B1 ) ( KSTAR1 ) ( CSTAR1 )  ' 

CALL  PRNT(KCC,NC12,NC12) 

CALL  UBOAT(Al, KCC, NC12, NC12, KOB1 ) 

PRINT*,'  A1  +  (B1 ) (KSTAR1 ) (CSTAR1 )  ' 

CALL  PRNT(KOBl, NC12, NC12) 

CALL  MMUL(BSTAR2, K5TAR2, NC22, E2, P2, BCG) 

CALL  MMUL (BCG, CSTAR2, NC22, P2, NC22, KCC) 

PRINT*,'  ( BSTAR2 ) ( KSTAR2 ) ( CSTAR2 )  ' 

CALL  PRNT ( KCC , NC22 , NC22 ) 

CALL  UBOAT(A2, KCC, NC22, NC22, K0B2) 

PRINT*,'  A2  +  (BSTAR2) (KSTAR2) (CSTAR2)  ' 

CALL  PRNT ( KOB2 , NC22 , NC22 ) 

IF( DEC. EQ. 3)  THEN 

CALL  MMUL(BSTAR3, KSTAR3,  NC32,  E3, NSEN, BCG) 

CALL  MMUL ( BCG, C3, NC32, NSEN, NC32, KCC ) 

PRINT*,'  ( BSTAR3 ) ( KSTAR3 ) ( C3 )  ’ 

CALL  PRNT ( KCC , NC32 , NC32 ) 

CALL  UBOAT(A3, KCC, NC32, NC32, K0B3) 

PRINT*,'  A3  +  (BSTAR3) (KSTAR3) (C3)  ' 

CALL  PRNT(K0B3, NC32,  NC32) 

ELSE 

CALL  MMUL(8STAR3, KSTAR3,  NC32,  E3, P3, BCG) 

CALL  MMUL (BCG, CSTAR3, NC32, P3,  NC32, KCC ) 

PRINT*,'  ( BSTAR3 ) ( KSTAR3 ) ( CST AR3 )  ' 

CALL  PRNT ( KCC , NC32 , NC32 ) 

CALL  UBOAT ( A3 , KCC , NC32 , NC32 ,  K0B3 ) 

PRINT*,'  A3  +  (BSTAR3) (KSTAR3) (CSTAR3)  ’ 

CALL  PRNT(KOB3, NC32,  NC32) 

CALL  MMUL(BSTAR4, KSTAR4,  NR2, E4, NSEN, BCG) 

CALL  MMUL ( BCG, C4, NR2, NSEN, NR2, KCC ) 

PRINT*,'  (BSTAR4) (KSTAR4) (C4)  ' 
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c* 


CALL  PRNT (KCC, NR2,  NR2) 

CALL  UB0AT(A4, KCC, NR2, NR2,  K0B4) 

PRINT*,1  A4  +  (B5TAR4) (KSTAR4) (C4)  1 

CALL  PRNT(K0B4,NR2,NR2) 

END  IF 

£5t!-tEI6RF(KOB1’  NC12>  NC0L>  W1- TEN,  NCOL,  STOR,  IER) 
PRINT*,'  EIGENVALUES  A1  +  (Bl) (KSTAR1 ) (CSTAR1 )  ' 
PRINT' (/) ' 

DO  421  I=1,NC12 

421  PRINT*,'  ' , W1 ( I ) 

PRINT' (/)' 

PRINT*,'  EIGENVALUES  A2  +  (BSTAR2) (K5TAR2) (CSTAR2) 
PRINT1 EIGRF(K0B2, NC22» NC0L’ 01 W1«  TEN, NCOL, STOR, IER) 

DO  422  1=1, NC22 

422  PRINT*,'  *  f  W1 (I ) 

PRINT' (/)' 

IF(DEC. EQ. 3)  THEN 

PRINT*,'  EIGENVALUES  A3  +  (BSTAR3) (KSTAR3) (C3)  ' 

ENDIF*' '  EIQENVALUES  A3  +  (BSTAR3)(KSTAR3)(CSTAR3) 

CALL  EIGRF (K0B3, NC32, NCOL, 8, Wl, TEN, NCOL, STOR, IER) 
PRINT’ (/)' 

DO  423  1=1, NC32 

423  PRINT*,'  ' , W1(I) 

PRINT' (/) ' 

IF(DEC. EQ. 4)  THEN 

PRINT*,'  EIGENVALUES  A4  +  (BSTAR4)(KSTAR4)(C4)  1 
CALL  EIGRF ( K0B4,  NR2,  NCOL , 8,  Wl ,  TEN, NCOL, STOR,  IER ) 
PRINT' (/) ' 

DO  424  1=1, NR2 

424  PRINT*,'  \W1(I) 

PRINT' (/)' 

ENDIF 


CALL  ABGC(A1, Bl, Cl, NC12, TKG1, NACT, NSEN, ABG1 ) 
PRINT*,'  A1  +  ( Bl ) ( TKG1 ) ( Cl )  ’ 

CALL  PRNT (ABG1, NC12, NC12) 

CALL  ABGC(A2, 82, C2, NC22, TKG2, NACT, NSEN, ABG2) 
PRINT*,'  A2  +  (B2)(TKG2)(C2)  ' 

CALL  PRNT( ABG2, NC22, NC22 ) 

CALL  ABGC(A3, B3, C3, NC32, TKG3,  NACT, NSEN, ABG3) 
PRINT*, '  A3  +  (B3) (TKQ3) (C3)  ' 

CALL  PRNT(ABG3, NC32, NC32) 

PRINT*,'  EIGENVALUES  OF  A1  +  (Bl ) (TKQ1 ) (Cl )  ' 

CALL  EIGRF ( A8G1, NC12, NCOL, 8, Wl, TEN, NCOL, STOR, IER) 
PRINT' (/)' 

DO  723  1=1, NC12 
723  PRINT*,'  \W1(I) 

PRINT' (/)' 

PRINT*,'  EIGENVALUES  OF  A2  +  ( B2 ) ( TKG2 ) ( C2 )  ’ 
CALL  EIGRF ( ABG2, NC22, NCOL, 8, Wl, TEN, NCOL, STOR, IER) 
PRINT' ( / ) ' 
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o  <r> 


DO  724  1=1, NC22 

724  PRINT*, '  ‘ . W1 ( I ) 

PRINT' (/)' 

PRINT*,'  EIGENVALUES  OF  A3  +  (B3) (TKG3) (C3)  ' 

CALL  £IGRF(ABG3, NC32, NCOL, 0, Wl, TEN, NCOL, STOR, IER ) 
PRINT1 (/)' 

DO  725  1=1, NC32 

725  PRINT*,'  \W1(I) 

PRINT' (/)' 

IF(DEC.EQ. 4)  THEN 

CALL  ABGC ( A4, 84, C4, NR2 , TKG4 , NACT , NSEN, ABG4) 
PRINT*,'  A4  +  (B4)(TKG4)(C4)  ' 

CALL  PRNT(ABG4,NR2,NR2) 

PRINT*,'  EIGENVALUES  A4  +  (B4)(TKG4)(C4)  ' 

CALL  EIGRF(ABG4,  NR2,  NCOL,  0,  Wl,  TEN,  NCOL,  STOR,  IER) 
PRINT’ (/)' 

DO  726  1=1, NR2 

726  PRINT*,'  ' , Wl ( 1 ) 

ENDIF 

PRINT' (/)' 

IF(NR. EQ. 0. OR. DEC. EQ. 4)  THEN 

GOTO  600 

ELSE 

DO  348  1=1, NACT 
DO  348  J=1,NS£N 
348  SAT(I,O)=0.0 

KT  =  -TKG 

CALL  UBOAT(SAT, TKG1,  NACT, NSEN, KT1 ) 

CALL  UBOAT ( SAT,  TKG2, NACT,  NSEN, KT2) 

CALL  UBOAT (SAT, TKG3, NACT, NSEN, KT3) 

DO  350  1=1, M 
DO  350  0=1, M 

350  MAJM(I, J)=0.0 

CALL  F0RMA(A1,D, W, NC1,NC12,  IC1) 

CALL  F0RMB(B1 , PH1A, NCI ,  NC12,  NACT, IC1 ) 

CALL  FORMCiCl, PHIS, NCI, NC12, NSEN, IC1 ) 

CALL  ABGC ( A1 , B1 ,  Cl,  NC12,  TKG1, NACT, NSEN, ABG5) 
PRINT*,'  A1  +  B1 ( TKG1 ) Cl  ' 

CALL  PRNT(ABG5,  NC12,  NC12) 

DO  351  1=1, NC12 
DO  351  J=1,NC12 

351  MA0M(I,J)=ABG5{I,0) 

CALL  F0RMA(A2,D,W,  NC2,NC22,  IC2) 

CALL  F0RMB(B2, PHI A, NC2, NC22, NACT, IC2) 

CALL  F0RMC(C2, PHIS,  NC2,  NC22,  NSEN, IC2) 

CALL  ABGC ( A2 , B2 , C2 , NC22 , TKG2 , NACT , NSEN , ABG6 ) 
PRINT*,'  A2  +  B2 ( TKG2 ) C2  ' 

CALL  PRNT(ABG6,NC22,NC  ’2) 

DO  352  I=1,NC22 
DO  352  J=l, NC22 

352  MAJM( I+K, 0+K)=ABG6(I, J) 
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1  0. 


call  FORMA (A3, 0,  W,  NC3,  NC32,  IC3) 

Call  F0RMB(B3,  FHIA,  NC3,  NC32.  r-JACT,  I C 3 ) 

CALL  F0RMC(C3, PHIS,  NC3,  NC32,  NSEN,  IC3) 

oStmtAB9° ( B3-''3VNC32’ TKG3<  NACT-  NSEN,  AB(j7 ) 
PRINT*,  A3  +  B3(TKG3)C3  ' 

CALL  PRNT(ABG7, NC32, NC32 ) 

DO  353  1=1, NC32 

DO  353  3=1, NC32 

353  MA3M(I+L, 3+L)=ABG7(I, J) 

CALL  MMUL(B2, KT1, NC22, NACT, NSEN, BCG) 

CALL  MMUL ( BCG, Cl ,  NC22,  NSEN,  NC12,  KCC ) 

PRINT*,'  B2(TKG1)C1  * 

CALL  PRNT(KCC, NC22, NC12) 

DO  354  1=1, NC22 
DO  354  J=l, NC12 

354  MA3M( I+K,  J)=KCC( I,  J) 

CALL  MMUL (B3.KT1, NC32,  NACT,  NSEN,  BCG) 

CALL  MMUL (BCG, Cl, NC32, NSEN,  NC12,  KOBl ) 

CALL  MMUL(B3, KT2, NC32, NACT, NSEN, BCG) 

CALL  MMUL (BCG, Cl, NC32, NSEN, NC12, K0B2 ) 

CALL  MADD(K0B1, K0B2, NC32, NC12, KCC) 

PRINT*,'  B3(TKG1)C1  +  B3(TKG2)C1  ' 

CALL  PRNT(KCC, NC32, NC12) 

DO  355  1=1, NC32 
DO  355  J=l, NC12 

355  MAJM( 1+L, J)=KCC( I, J) 

CALL  MMUL(B3, KT2, NC32, NACT,  NSEN,  BCG) 

CALL  MMUL (BCG, C2, NC32,  NSEN,  NC22,  KCC) 

PRINT*, '  B3(TKG2)C2  ' 

CALL  PRNT(KCC, NC32,  NC22) 

DO  356  1=1, NC32 
DO  356  3=1, NC22 

356  MAJM( I+L, 3+K)=KCC(1, 3) 

CALL  MMUL (B4, KT1,  NR2,  NACT,  NSEN, BCG) 

CALL  MMUL (BCG, Cl, NR2, NSEN, NC12, KOBl ) 

CALL  MMUL (B4, KT2, NR2,  NACT,  NSEN,  BCG ) 

CALL  MMUL ( BCG, Cl , NR2 , NSEN, NC12 , K0B2 ) 

CALL  MMUL(B4, KT3,  NR2,  NACT,  NSEN,  BCG) 

CALL  MMUL ( BCG, C 1 , NR2 ,  NSEN , NC 1 2 ,  K0B3 ) 

CALL  MADD(K0B1, K0B2,  NR2,  NC12,  BCG) 

CALL  MADD ( BCG , K0B3 , NR 2 , NC 1 2 , K0B4 ) 

PRINT*,'  B4(TKG1)C1  +  B4(TKG2)C1  +  B4(TKG3)C1 
CALL  PRNT ( K0B4 , NR2 , NC12 ) 

DO  357  1=1, NR2 
DO  357  3=1,  NC12 

357  MAJM(I+P5,3)=K0B4(I,3) 

CALL  MMUL(B4, KT2, NR2, NACT, NSEN, BCG) 

CALL  MMUL (BCG, C2, NR2, NSEN, NC22, K0B2) 

CALL  MMUL (B4,KT3,NR2, NACT, NSEN, BCG) 

CALL  MMUL (BCG, C2, NR 2, NSEN, NC22, K0B3) 

CALL  MADD(K0B2, K0B3, NR2, NC22, K0B4) 

PRINT*,'  B4(TKG2)C2  +  B4(TKG3)C2  ' 

CALL  PRNT(K0B4, NR2, NC22) 
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58 

DO  858  0=1,  NC22 

358  MAOM(  I+P5,  J+K )  =K0B4(  1 ,  ,1) 

CALL  MMUL(B4, K] 3, NR 2, NACT, N5EN, KOBl 1 
CALL  MMUL  ( KOBl ,  C'3,  NR 2,  N5EN,  NC32,K0B2 ) 

PRINT*,’  B4(TKG3)C3  1 
CALL  PRNT ( K0B2, NR2, NC32 ) 

DO  359  1=1, NR2 
DO  359  J=l, NC32 

359  MAJM(I+P5, J+L)=K0B2(I,  J) 

CALL  MADD(KT1,KT2,  NACT,  NSEN,  BCG) 

CALL  MADD(BCG,KT3, NACT, NSEN,  KCC) 

CALL  MMUL ( KCC, C4, NACT, NSEN, NR2, BCG) 

CALL  MMUL(B1, BCG, NC12, NACT, NR2, ABK1 ) 

PRINT*, 1  B1 ( TKG1  +  TKG2  +  TKG3)C4  * 

CALL  PRNT (ABK1, NC12,  NR2) 

DO  361  1=1, NC12 
DO  361  J=l, NR2 

361  MAJM(I, J+P5)=ABK1(I,  J) 

CALL  MMUL (B2, BCG, NC22, NACT, NR2.ABK2) 

PRINT*,1  B2(TKG1  +■  TKG2  +  TKG3)C4  1 
CALL  PRNT (ABK2, NC22, NR2 ) 

DO  362  1=1, NC22 
DO  362  J=l, NR2 

362  MAJM(I+K,0+P5)=ABK2(I,0) 

CALL  MMUL ( B3 , BCG, NC32 ,  NACT , NR2 ,  ABK 3 ) 

PRINT*, 1  B3(TKG1  +  TKG2  +  TKG3)C4  1 
CALL  PRNT(ABK3, NC32,  NR2) 

DO  363  I=1,NC32 
DO  363  0=1, NR2 

363  MAJM(I+L, 0+P5)=ABK3( I, J) 

CALL  MADD(TKG1,  TKG2, NACT, NSEN,  BCG) 

CALL  MADD( BCG, TKG3, NACT,  NSEN,  KCC ) 

CALL  F0RMA(A4, D, W, NR, NR2, IR) 

CALL  F0RMB(B4, PHIA, NR, NR2, NACT, IR) 

CALL  F0RMC(C4, PHIS, NR,  NR2,  NSEN,  IR) 

CALL  ABGC(A4, B4, C4, NR2,  KCC,  NACT,  NSEN,  ABG8) 

PRINT*, 1  A4  +  B4( TKG1  +  TKG2  +  TKG3 )C4  ' 

CALL  PRNT (ABG8, NR2, NR2) 

DO  365  1=1, NR2 
DO  365  J=l, NR2 

365  MAJM(I+P5,>P5)=ABG8(I,0) 

PRINT*, 1  LOWER  TRIANGLE  THREE  CONTROLLER  MAJOR  MATRIX 
PRINT*,  ’  WITH  RESIDAIJLS 

CALL  PRNT(MAJM,M,  M) 

PRINT’ (///)' 

PRINT*, '  TRANSFORMED  SYSTEM  EIGENVALUES  ’ 

CALL  EIGRF(MAJM, M, NDA, 8, 2, TEN, NCOL,  WORK, IER) 

PRINT’ ( /) ' 

DO  366  1=1, M 

366  PRINT*,'  ■ , Z( I) 

PRINT' (//) ' 


c~>  <~j  o  o  o  <->  o 


GW  CONTINUE 


PRINT* (///)' 

PRINT*’ '  ** 

PRINT*, '  **  ENO  OOFB335  PROGRAM 
PRINT*, 1  ** 

PRINT*,  1  ************  ***************** 

PRINT' (///) ' 

ENO 

SUBROUTINE  FORMX2(XO,  INIT) 


FORMX2  IS  THE  SAME  FOR  THREE  OR  FOUR 
CONTROLLERS  SINCE  THE  RESIDUALS  BECOME 
THE  FOURTH  CONTROLLER. 

COMMON/NUM/IC1 ( 21 ) , IC2 ( 21 ) , IC3 ( 21 ) , IR ( 21 ) , NCI , NC2, NC3 , NR 
REAL  X0(44),INIT(4,12) 

INTEGER  M,I,J,K,L 
DO  1  1=1, NCI 
M  =  IC1(I) 

XO(I)  =  INIT(1,M) 

XO(I+NCl)  =  INIT(2,M) 

X0(I+NC1*2)  =  INIT(3,M) 

1  XO(I+NCl*3)  =  INIT(4,M) 

J  =  NC1*4 

DO  2  1=1, NC2 
M  =  IC2 ( I ) 

X0(I+0)  =  INIT(1,M) 

XO(I+J+NC2)  =  INIT  (2,  M) 

X0(I+J+NC2*2)  =  INIT ( 3, M) 

2  XO(I+J+NC2*3)  =  INIT(4,M) 

K  =  J  +  NC2*4 

DO  3  1=1, NC3 
M  =  IC3( I ) 

XO(I+K)  =  INIT ( 1, M) 

XO(I+K+NC3)  =  INIT(2,M) 

XO(I+K+NC3*2)  =  INIT ( 3, M) 

3  XO(I+K+NC3*3)  =  INIT (4,  M) 

L  =  K  +  NC3*4 

DO  4  1=1,  NR 
M  =  IR(I) 

XO(I+L)  =  INIT(1,M) 

XO(I+L+NR)  =  INIT ( 2, M) 

XO(I+L+NR*2)  =  INIT(3,M) 

4  XO(I+L+NR*3)  =  INIT(4,M) 

END 

SUBROUTINE  FACTOR (N, A, S, MR) 

C  A=S'S 

DIMENSION  A(1),S{1) 
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COMMON/MAINB/  NCOL.NCOLl 
COMMON/ I NOUT/KOUT 
T0L=l.E-6 
MR=0 

NN=N*NCOL 

TOL1=0. 

DO  1  I=l,NN,NCOLl 
R=ABS(A(I)) 

1  IF  (R.OT.TOL1)  TOLl=R 
T0L1=T0L1*1.E-12 
11=1 

DO  50  1=1, N 
IM1=I-1 

DO  5  JJ=I,NN,NCOL 
5  S(OO)=0. 

ID=II+IM1 

R=A( ID) -DOT ( IM1 ,S(II)fS(II)) 

IF  (ABS(R).LT.  (TOL*A(ID)+TOLl))  GO  TO  50 
IF  (R)  15,50,20 
15  MR=-1 

WRITE (KOUT, 1000) 

1000  FORMAT(37H0TRIED  TO  FACTOR  AN  INDEFINITE  MATRIX 
RETURN 

20  S(ID)=SQRT(R) 

MR=MR+1 

IF  (l.EQ.N)  RETURN 

L=II+NCOL 

DO  25  00=L,NN,NC0L 

I0=00+IM1 

25  5(I0)=(A(I0)-D0T(IM1,S(II),S(00)))/S(ID) 

50  II=Il+NCOL 
RETURN 
END 

SUBROUTINE  FORMA(A,D,W,N,N2,IC) 

COMMON/MAINB/NCOL 

REAL  A(NCOL,NCOL), W(17),D(17) 

INTEGER  IC(N), I, J, N,M 

DO  1  1=1, N2 

DO  1  0=1, N2 

A(I,J)=0.0 

CONTINUE 

DO  2  1=1, N 

M=  IC(I) 

A((I+N), (I+N))=D(M) 

A(I, (I+N) )  =  1.0 
A((I+N), I)  =  -(W(M)**2) 

!  CONTINUE 
RETURN 
END 

SUBROUTINE  FORMS (B, PHI , N, N2 , NACT, IC ) 

COMMON/MAINB/NCOL 

REAL  B(NCQL,NCOL),PHI(NCOL,NCOL) 

INTEGER  IC(N), NACT, N,M, 1,0 


DO  1  1=1,  N2 
DO  1  0=1,  NACT 
B(  I,  J)  =  0.0 
CONTINUE 
DO  2  1=1, N 
M  =  IC(I) 

DO  2  J=1,NACT 
B((N+I), J)  =  PHI (M, 0) 

CONTINUE 

RETURN 

END 

SUBROUTINE  FORMC(C, PHIS, N,N2,  NSEN,  IC) 

COMMON/MAINB/NCOL 

REAL  C(NCOL,NCOL), PHIS(NCOL, NCOL) 

INTEGER  IC(N) , M, NSEN, N, N2, I, J 

DO  1  1=1,  NSEN 

DO  1  0=1,  N2 

C(I,0)  =  0.0 

CONTINUE 

DO  2  1=1, NSEN 

DO  2  0=1, N 

M  =  IC(O) 

C(I,N+0)=PHIS(M, I) 

CONTINUE 

RETURN 

END 

SUBROUTINE  F0RMQ(Q,  A,N, IC) 

COMMON/MAINB/NCOL 

REAL  h( NCOL ),Q( NCOL, NCOL) 

INTEGER  1,0, K,  M, N, N2, IC(NCOL) 

N2  =  N  *  2 
DO  1  1=1, N2 
DO  1  0=1,  N2 
Q( 1 , 0)  =  0.0 
CONTINUE 
DO  2  1=1, N 
0  =  I 
M  =  IC(I) 

DO  3  K=I— 1, I 
Q(I+K, 0+K)  =  A(M) 

CONTINUE 

RETURN 

END 

SUBROUTINE  GMINV(NR,NC, A,U,MR,MT) 
DIMENSION  A(1),U(1) 

COMMON/MAIN1/  NDIM, NDIM1, S(1 ) 

COMMON/MAINB/NCOL, NCOL 1 

COMMON/ INOUT/KOUT 

T0L=1. E-12 

MR=NC 

NRM1=NR-1 

T0L1=1. £-20 


O' 


w 


DO  m  J=l, NC 
FAC=DOT(NR,A(OJ), A(JJ) ) 

OMl=J-l 
ORM=JJ+NRMl 
JCM=JJ+JM1 
DO  20  I=JJ, JCM 
20  U(I)=0. 

U(JCM)=1.0 

IF  (J.EQ.l)  60  TO  54 

KK=1 

DO  30  K=l, JM1 
IF  (S(K).EQ. 1.0)  GO  TO  30 
TEMP=— DOT ( NR , A ( JO ) , A ( KK ) ) 

CALL  VADD(K,TEMP,U( JJ),U(KK)) 

30  KK=KK+NCOL 
DO  50  L=l,  2 
KK=1 

DO  50  K=1,0M 
IF  (S(K).EQ.0. )  GO  TO  50 
TEMP=-DOT(NR, A( JO), A(KK) ) 

CALL  VADD(NR,TEMP, A( JO), A(KK) ) 

CALL  VADD(K,TEMP,U(JJ),U(KK)) 

50  KK=KK+NCOL 
T0L1=T0L*FAC 
FAC=DOT(NR,A(JJ),A(JO)) 

54  IF  (FAC. GT. T0L1)  GO  TO  70 
DO  55  I=JJ, JRM 

55  A(I)*0. 

S(J)=0. 

KK*1 

DO  65  K=l, JM1 

IF  (S(K). EQ.0. )  GO  TO  65 

TEMP=-DOT(K, U(KK) ,  U( JJ) ) 

CALL  VADD(NR,TEMP,A(JO), A(KK)) 

65  KK=KK+NCOL 

FAC=DOT(J,U(OJ),U(JO)) 

MR=MR-1 
GO  TO  75 
70  S(J)=1.0 

KK=1 

DO  72  K=1,0M1 

IF  (S(K).EQ.l.)  GO  TO  72 

TEMP=-DOT(NR,A(JO),A(KK)) 

CALL  VADD(K,TEMP,U(JJ),U(KK)) 

72  KK=KK+NCOL 
75  FAC* 1./ SQRT (FAC) 

DO  80  1=00, JRM 
80  A(I)*A(I)*FAC 
DO  85  I=JJ, JCM 
85  U(I)*U(I)*FAC 
100  JO=JJ+NCOL 

IF  (MR.EQ.NR.OR.MR.EQ.NC)  GO  TO  120 
IF  (MT.NE.0)  WRITE(KOUT, 110)NR, NC,MR 
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110  FORMAT ( 13, 1HX, 12, 8H  M  RANK, 12) 

126  NENO=NONCOL 
30=1 

DO  135  0=1, NC 
DO  125  1=1, NR 
II=I-J 
S(I)=0. 

DO  125  KK=00,  NEND, NCOL 
125  S(I)=S(I)+A(II+KK)*U(KK) 

11=0 

DO  130  1=1, NR 
U(II)=S(I) 

130  II=II+NC0L 
135  00=J0+NC0L1 
RETURN 
END 

SUBROUTINE  GMINV1(NR,NC,A,U,MR,MT) 
DIMENSION  A(1),U(1) 

C0MM0N/MAIN1/  NDIM, NDIM1, S{ 1) 

COMMON/MAINB/NCOL,  NC0L1 

COMMON/MAINA/NDA, NDA1 

COMMON/ INOUT/KOUT 

T0L=1. E-12 

MR=NC 

NRM1=NR-1 

T0L1=1. E-20 

0J=1 

DO  100  0=1, NC 
FAC=DOT(NR,A(00),A(00)) 

0M1=0-1 
0RM=00+NRM1 
0CM=0O+0Ml 
DO  20  1=00, OCM 
20  U(I)=0. 

U( OCM) =1.0 

IF  (O.EQ.l)  GO  TO  54 
KK=1 

DO  30  K=l, 0M1 
IF  (S(K).EQ.1.0)  60  TO  30 
TEMP=-DOT(NR,A(00),A(KK) ) 

CALL  VADD(K, TEMP, U(00), U(KK) ) 

30  KK=KK+NDA 
DO  50  L=l, 2 

KK=1 

DO  50  K=l, 0M1 

IF  (S(K).EQ.e. )  GO  TO  50 

TEMP=-D0T(NR,A(00),A(KK)) 

CALL  VADD(NR, TEMP, A(00), A(KK) ) 
CALL  VADD(K, TEMP, U(00), U(KK) ) 

50  KK=KK+NDA 
T0L1=T0L*FAC 
FAC=D0T(NR,A(00),A(00)) 

54  IF  (FAC.GT.T0L1)  60  TO  70 
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DO  55  I=JJ, JRM 
55  A( I)=0. 

S(J)=0. 

KK=1 

DO  65  K=l, JM1 

IF  (S(K).EQ.e.)  GO  TO  65 

TEMP=-DOT(K,U(KK),U(JJ)) 

CALL  VADD ( NR , TEW, A{ JJ), A(KK)) 

65  KK=KK+NDA 

FAC=DOT ( J,  U(  JJ) ,  U(  JJ ) ) 

MR=MR-1 
GO  TO  75 
70  S(J)=1.0 

KK=1 

DO  72  K=l, JM1 
IF  (S(K).EQ.l.)  GO  TO  72 
TEMP=-DOT(NR, A( JJ), A(KK)) 

CALL  VADD(K, TEMP,U( JO),U(KK)) 

72  KK=KK+NDA 
75  FAC=1. /SQRT(FAC) 

DO  80  I=JJ, JRM 
80  A(I)=A(I)*FAC 
DO  85  I=JJ, JCM 
85  U(I)=U(I)*FAC 
100  JJ=JJ+NDA 

IF  (MR.EQ.NR.OR.MR.EQ.NC)  GO  TO  120 
IF  (MT.NE.0)  WRITE(KOUT,110)NR,NC,MR 
110  FORMAT ( 13, 1HX, 12, 8H  M  RANK, 12) 

120  NEND=NC*NDA 
JJ=1 

DO  135  0=1, NC 
DO  125  1=1, NR 
II=I-J 
S(I)=0. 

DO  125  KK=00,  NENO,  NDA 
125  S(I)=S(I)+A(II+KK)*U(KK) 

II=J 

DO  130  1=1, NR 
U(II)=S(I) 

130  II=II+NDA 
135  J0=00+N0A1 
RETURN 
END 

SUBROUTINE  INTEG(N,A,C,S,T) 

C  S=INTEGRAL  EA*C*EA  FROM  0  TO  T 

C  C  IS  DESTROYED 

DIMENSION  A(1),C(1),S(1) 
COMMON/MAIN1/  NDIM, NDIM1,  X(l) 
COMMON/MAINB/NCOL, NC0L1 
COMMON/MAIN2/COEF ( 100 ) 

NN=N*NCOL 

NM1=N-1 

IND=0 
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ANORM=  XNORM ( N ,  A ) 

DT=T 

5  IF  (ANORM*ABS(DT).LE.0.5)  GO  TO  10 
DT=DT/2. 

IND=IND+1 
GO  TO  5 

10  00  15  1=1,  NN,  NCOL 

0=I+NM1 
DO  15  00=1,0 
15  S(  00)  =DT*C(  30) 

Tl=DT+*2/2. 

00  25  IT=3, 15 

CALL  MMUL(A, C,N, N,  N,  X) 

00  20  1=1, N 
II=(I-1)*NC0L 
DO  20  00=1, NN, NCOL 
11=11+1 

C(00)=(X(00)+X(II))*T1 
20  S(30)=S(00)+C(00) 

25  T1=0T/FL0AT(IT) 

IF  (INO.EQ.0)  GO  TO  100 
COEF( 11 )=1 . 0 
DO  30  1=1,10 
11=11-1 

30  COEF(II) =DT*COEF ( I 1+1 ) / FLOAT ( I ) 

11=1 

fT*  DO  40  1=1,  NN,  NCOL 

0=I+I\M1 
DO  35  00=1,0 
35  X(00)=A(00)*C0EF(1) 

X(II)=X(II)+C0EF(2) 

40  I1=II+NC0L1 
DO  55  L=3, 11 
CALL  MMUL(A, X,N, N, N, C) 

11=1 

T1=C0EF(L) 

DO  55  1=1, NN, NCOL 
0=I+NM1 
DO  50  JJ=I,J 
50  X(00)=C(00) 

X(  II )=X( II)+T1 
55  II=II+NC0L1 
C  X=EXP(A*DT) 

L=0 

60  L=L+1 

CALL  MMUL(X,S, N,N, N, C) 

11=1 

DO  90  1=1, N 
0=11 

IF  (I.EQ.l)  GO  TO  75 
DO  70  30=1, II, NCOL 
S(OJ)=S(J) 

70  0=3+1 
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I 


I 

I 

I 


•  J 


80 

85 

87 

90 


100 


5 


DO  85  00=1, N 
KK=OJ 

DO  80  K=I,NN, NCOL 
S(0)=S(0)+C(K)*X(KK) 
KK=KK+NCOL 
0=0+NC0L 

DO  87  00=1, NN, NCOL 
C(00)=X(JJ) 

II=II+NCOL 

IF  (L.EQ.IND)  GO  TO  100 

CALL  MMUL(C,C,N,  N,  N,  X) 

GO  TO  60 

CONTINUE 

RETURN 

ENO 


SUBROUTINE  MEXP(N, A, T, EA) 

DIMENSION  A( 3364) , EA(3364) , C( 100),  D( 101 ) , E( 100) 
COMMON/MAINA/NDA, NDA1 


C0MM0N/MAIN1/N0IM, ND1M1, TEN, X(3364) 
COMMON/MAINB/NCOL , NC0L1 
NN=NDA*N 


NM1=N-1 

PRINT*, 1  MEXPIN  ' 

IF  (N.GT.l)  GO  TO  5 

EA(1)=EXP(T*A(1)) 

RETURN 


W=1.0 


DO  10  1=1, NN, NDA 
1L=I+NM1 
DO  10  J=I, IL 
10  EA( J)=A( J) 
C1=XN0RM1(N, A) 
IND=0 
L=1 


T1=T 


15  IF  (ABS(T1*C1).LE.3.0)  GO  TO  20 
Tl=Tl/2. 


IND=IND+1 
GO  TO  15 
20  C2=0. 


DO  25  1=1, NN, NDA1 
25  C2=C2-EA(I) 
C2=C2/FL0AT(L) 
C(L)=C2 
D(L+1)=0. 

II=N+1-L 
E( II )=W 
11=1 


DO  35  1=1, NN, NDA 
IL=I+NM1 
DO  30  J=I, IL 
30  X(0)=EA(0) 

X( II)=X(II)+C2 
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35  I I =1 I+NDA1 

IF  (L.EQ.N)  GO  TO  40 
CALL  MMUL(X,A,N,  N,  N,  EA) 
W=W*T1/FL0AT(L) 

L=L+1 
GO  TO  20 
40  CONTINUE 

C ******  CAN  CHECK  X  0  FOR  ACCURACY 
0=N+25 
DO  50  L=N,  J 
DO  45  K=l, N 

D(K)=(D(K+1)-W*C(K) ) *T 1 / FLOAT ( L ) 

45  E(K)=E(K)+D(K) 

50  W=D(1) 

11=1 

DO  60  1=1, NN, NDA 
IL=I+NM1 
DO  55  J=I, IL 
55  EA(0)=E(1)*A(0) 

EA(1I)=EA(II)+E(2) 

60  II=II+NDA1 

IF  (N.EQ.2)  GO  TO  85 
DO  80  L=3, N 

CALL  MMUL(EA,  A,N,N,N, X) 

11=1 

C2=E(L) 

DO  75  1=1,  NN,  NDA 
IL=I+NM1 
DO  70  0=1, IL 
70  EA(J)=X(J) 

EA(II)=EA(II)+C2 
75  II=II+NDA1 
80  CONTINUE 

85  IF  (IND.EQ.0)  RETURN 
DO  100  L=l, IND 
DO  90  1=1,  NN,  NDA 
IL=I+NM1 
DO  90  J=I, IL 
90  X(0)=EA(0) 

100  CALL  MMUL ( X, X, N, N, N, EA) 

RETURN 

END 

SUBROUTINE  MLINEQ(N,  A,  C, X, TOL,  IER) 
C  SOLVES  A1 X+XA+C=0 

C  A  AND  X  CAN  BE  IN  SAME  LOCATION 
C  ANSWER  RETURNED  IN  C  AND  X 

DIMENSION  A(1),C(1),X(1) 
COMMON/MAINB/  NCOL,  NC0L1 
C0MM0N/MAIN3/F(1) 

ADV=T0L*1.  E-6 
DT=.  5 
DT1=0. 

NN=N*NC0L 


5 


15 

20 


C 


25 

30 

35 

40 


60 


65 

70 


75 


80 


DO  5  11=1, NN.NCOLl 
0T1=DT1-A(1I) 

DU=DT1/N 

IF  (DTI.  GT. 4.0)  DT=DT*4. 0/DT1 
11=1 

DO  20  1=1, N 
DO  15  00=1, NN, NCOL 
X(00)=DT*A(00) 

X( II )=X( II ) -. 5 
II=II+NC0L1 

CALL  GMINV(N, N,X,F,MR,0) 

IER=4 

IF  (MR.NE.N)  RETURN 
CALL  MMUL(C,F,N,N,  N,  X) 

INITIALIZATION  OF  X,  F 
1=1 

DO  40  11=1, NN, NCOL 
0=11 

IF  (I.EQ. 1)  GO  TO  30 
DO  25  00=1, II, NCOL 
C(0)=C(00) 

0=0+1 

ID=0 

DO  35  00=11, NN, NCOL 
C(0)=DT*D0T(N,F(II),X(00)) 

0=0+1 

F(ID)=F( ID)+1. 0 
1=1+1 

DO  90  1T=1, 20 
NEZ=0 

CALL  MMUL(C, F, N, N, N, X) 

11=1 

0=1 

GO  TO  70 
0=11 

DO  65  00=1, II, NCOL 
C(0)=C(00) 

0=0+1 

ID=0 

DT1=C(0) 

DO  75  00=II,NN, NCOL 
C(0)=C(0)+D0T(N, F(II), X(00) ) 

0=0+1 

0=0-1 

DO  80  00=11,0 
X(00)=F(00) 

IF  (ABS(C(ID) ). GT. 1. E150 )  GO  TO  95 
I Fr ( ABS ( C ( ID ) -DTI ) . LT. ( ADV+TOL*ABS (C( ID) ) ) )  NEZ=NEZ+1 

II=II+NCOL 

IF  (I.LE.N)  GO  TO  60 
IF  (NEZ.EQ.N)  GO  TO  150 
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CALL  MMUL(X,X,N,N,N,  F) 

90  CONTINUE 
95  IER=1 
RETURN 

150  CONTINUE 
NM1=N-1 

DO  155  1=1, NN, NCOL 
II=I+NM1 
DO  155  JJ=I,II 
155  X(JJ)=C(JJ) 

IER=0 

RETURN 

END 

SUBROUTINE  MMUL(X,  Y,N1,N2,N3,Z) 

COMMON  /MAINB/  NCOL 

DIMENSION  X(NCOL,  1 ), Y(NCOL,  1 ),  Z{NCOL, 1 ) 
DO  3  J=1,N3 
DO  2  1=1,  N1 
S=0. 

DO  1  K=1,N2 

1  S=S+X(I,K)*Y(K,J) 

2  Z(I,0)=S 

3  CONTINUE 
END 

SUBROUTINE  MRIC(N, A, S, Q, X, Z, TOL, IER) 
DIMENSION  A(1),S(1),Q(1),X(1),Z(1) 
COMMON/MAIN1/NDIM,  NDIM1,  F(l) 
COMMON/MAINB/NCOL, NCOL1 
C0MM0N/MAIN2/TR( 1 ) 

COMMON/ INOUT/KOUT 

ADV=TOL*l.  E-6 

NN=N*NCOL 

NM1=N-1 

IND=1 

COUNT=0. 

IF  (IER. EQ. 1)  C0UNT=99. 

IF  (IER. EQ. 1)  MR=N 
IF  (IER. EQ. 1)  GO  TO  100 
Tl=-1. 

300  CONTINUE 
IER=0 

COUNT=COUNT+l. 

DO  15  1=1, N 
DO  15  J=I,  NN,  NCOL 
15  X( J)=-S(0) 

CALL  INTEG(N, A, X, Z, T1 ) 

CALL  FACTOR (N,Z,X, MR) 

IER=1 

IF  (MR.LT.0)  GO  TO  200 
IER=0 

CALL  GMINV(N, N, X,Z,MR,0) 

CALL  TFR(TR,Z,N,N, 1,2) 

CALL  MMUL(Z, TR, N, N,  N,  X) 


*  T'«~r  -  r 


1 


'J* 


I 


El  frl 


M 


I# 


I .  - 

L-  * 

i  *  .. 


DO  18  11=1, NN, NCCL1 
1  =  11 


;ci  "■ 

DO  17  0=11 , NN,  NCOL 
X(0)=(X(0)+A(I))/2 

X(1)=X(0) 

• 

17 

1=1+1 

K  ' 

• 

18 

CONTINUE 

r-  . 

100 

CONTINUE 

fa 

DO  16  1=1, N 

r 

16 

TR(I)=-1.0 

C 

A+SX  IS  STABLE 

20 

101 


25 


30 

35 

40 

45 

50 

55 


POSSIBLE  UNCONTROLLABILITY  IF  MR.NE.N 
TIM  DILLQW  IS  A  NUTTY  MATH  PROF 
TOL1=TOL/10. 

MAXIT=40 

DO  40  IT=1,MAXIT 
IF  (IER.EQ.l)  GO  TO  101 
CALL  MMUL(S,X,N,N,N,F) 

CALL  MMUL(X,F,N,N,N,Z) 

DO  20  1=1, NN, NCOL 
II=I+NM1 
DO  20  0=1,  II 
X(0)=A(0)-F(0) 

Z( J)=Z( J)+Q(U) 

CONTINUE 
IER=0 

CALL  MLINEQ(N, X,  Z,  X,  T0L1,  IER) 

IF  (IER.NE.0)  GO  TO  200 
L=0 
C1=0.  0 
11  =  1 

DO  25  1=1, N 

IF  (ABS(X(II)-TR(I) ).LT. (ADV+TOL*X( II ) ) 

TR( 1 )=X( II ) 

II=II+NC0L1 
C1=C1+TR( I j 

IF  (ABS(Cl).GT. 1. E20)  GO  TO  50 
IF  (L.NE.N)  GO  TO  40 
CALL  GMINV(N1N,Z,F,MR,0) 

CALL  MMIJL(S,X,N,N,N,Z) 

DO  30  1=1, NN, NCOL 
II=I+NM1 
DO  30  J=I,II 
Z(0)=A(0)-Z( J) 

IF  (MR.NE.N)  WRITE(K0UT,35)MR 

FORMAT (27H0RICCATI  SOLN  IS  PSD — RANK  ,13) 

GO  TO  65 
CONTINUE 

WR I T E ( K OUT , 45)  MAX1T 

FORMAT ( 27H0RICCATI  NON-CONVER'.jENT  IN  ,I2,11H  ITERATIONS) 
GO  TO  60 

W  OE(KOUT,55)IT,T1 

FORMAT(30H0RICCATI  BLOW-UP  AT  ITERATION  ,I2,12H  INITIAL  T-- 
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L=L+1 


, F10. 5) 


OH  1ER=1 
o5  RETURN 

2m  IF  (IND.EQ.2)  GO  TO  250 
IF  (COUNT. GE. 10. )  RETURN 
T1=T1/(2.**C0UNT) 

IND=2 
GC  TO  300 

250  T1=T1*(2.**C0UNT) 

IND=1 
GO  TO  300 
END 

SUBROUTINE  PRNT(MAT,N,M) 
COMMON/MAINB/NCOL 
REAL  MAT ( NCOL , NCOL ) 

INTEGER  N,  I,  J,  K, M 
PRINT*. 1  1 
IF  (M.GT. 12)  GOTO  2 
DO  1  1=1, N 

PRINT' ( IX, 12F10. 4) 1 , (MAT ( I, J), J=1,M) 

1  CONTINUE 
GOTO  10 

2  CONTINUE 

IF  (M.GT. 24)  THEN 
CALL  PRNTXL(MAT,N,M) 

RETURN 
END  IF 

DO  3  1=1, N 

PRINT'  ( IX,  12F10. 4) ' , (MAT ( I, J), J=l, 12) 

3  CONTINUE 
PRINT' (//)' 

DO  4  1=1, N 

PRINT' (IX, 12F10, 4) ' , (MAT(I, J ) , J= 1 3 , M ) 

4  CONTINUE 

10  PRINT '(///)' 

RETURN 

END 

SUBROUTINE  PRNTXL(MAT, N,M) 
COMMON/MAINB/NCOL 
COMMON/MAINA/NDA 
REAL  MAT(NDA.NDA) 

INTEGER  I, J, K, L,M, N 
PRINT*, '  ' 

DO  1  L=l, M, 12 
K  =  L  +  11 

IF  (M-L.LT.ll)  K  =  M 
DO  2  1=1, N 

PRINT ' ( IX, 12F10. 5) ', (MAT ( I , J) , J=L, K) 

2  CONTINUE 

PRINT '(//)' 

1  CONTINUE 

PRINT’ (///) ’ 

RETURN 

END 


1 


i 

i 

i 


X 

4 

1 
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SUBROUTINE  TFR; X, A, N, M, K, I ; 


GIVES 

GIVES 

GIVES 


X  =  A 
X  =  A 
X  =  A 


A  AS  A  VECTOR 


4  GIVES  A  =  X  WHERE  X  WAS  A  VECTOR 

DIMENSION  X(1),A(1) 

COMMON/MAINB/NCOL 

JS=(K-1 )*NCOL*M 

JEND=M*NCOL 

GO  TO  (10, 30, 50, 70), I 

DO  20  11=1, N 

DO  20  -3-3=11, JEND.NCOL 

X(0J)=A(00+0S) 

RETURN 

DO  40  11=1, N 
KK=(1I-1)*NC0L 
DO  40  -3-3=1,  M 
LL=( JJ-1 )*NCOL+II 
X(KK+-3J)=A(LL+JS) 

RETURN 

KK=0 

DO  60  11=1, OEND, NCOL 
LL-'II+N-l 
DO  60  0-3=11,  LL 
KK=KK+1 

X  ( KK )  =A  ( J-3+JS ) 

RETURN 
KK=M*M+1 
DO  80  11=1,  M 
LL=(M-II)*NC0L+1 
DO  80  1-3=1,  N 
KK=KK-1 
J-3=LL+N-IJ 
A(  -3-3+-3S)=X(KK) 

RETURN 

END 

SUBROUTINE  TIME ( EAT2, MM, DT, XI, XO, MODE, EAT, WORK, DEC ) 
C0MM0N/MA1NA/NDA, NDAl 
COMMON/SAVE/T ( 100 ) , TS( 100 ) 

C0MM0N/Nlfi/ICl(21 ), IC2 ( 21 ) , IC3 ( 21 ) , IR(21 ), NCI, NC2, NC3, NR 

REAL  XO ( NDA ) , EAT2 ( NDA , NDA ) , EAT ( NDA , NDA ) , DT , MODE ( 2 , 1 2 ) 

R EAL  WORK ( NDA , NDA ) , A , AA , Z , X 1 ( ND A ) 

INTEGER  DEC,  I,  J,K,L,M,N,KK,LL,MM 

PRINT*, 1  TIMEIN  1 

00  222  1=1, MM 

DO  222  -3=1,  MM 

WORK  (I,  -3)  =  EAT2(  I,  -3) 

N  =  1 
Z  =  DT*5. 

PRINT*, 1  TIME  X  Y  1 

KK  =  0 
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A  =  0 

2D3  CONTINUE 
M  =  0 

201  CONTINUE 
PRINT*,  *  201  1 

IF  (KK.GT.0)  THEN 

CALL  VMULFF(WORK,  EAT2,  MM,  MM,  MM,  NDA,  NDA,  EAT,  NDA,  IER ) 
PRINT*, 1  201  V  1 
DO  112  1=1, MM 
DO  112  J=1,MM 
EAT2(I,0)  =  EAT ( I, J) 

112  CONTINUE 
PRINT*,'  112  ' 

ENDIF 

KK  =  1 

CALL  VMULFF(EAT2, XO,MM,MM, 1, NDA, NDA, XI, NDA, IER) 
PRINT*, 1  112  V  ' 

DO  113  1=1, MM 

113  XO(I)  =  X1(I) 

M  =  M+l 

IF  ( (M*DT). EQ. Z)  THEN 
A  =  A  +  2 
DO  202  K=l, 2 
AA  =  0.0 
DO  204  1=1, NCI 
0  =  IC1(I) 

204  AA  =  AA  +  MODE{K,J)  *  XI ( I ) 

LL  =  NCI  *  4 

DO  205  1=1, NC2 
J  =  IC2(I) 

205  AA  =  AA  +  MODE(K,J)  *  X1(I+LL) 

L  =  KK  +  NC2  *  4 

DO  206  1=1, NC3 
J  =  IC3(I) 

206  AA  =  AA  +  MODE(K, J)  *  X1(I+L) 

IF  (DEC.EQ.4)  THEN 

M  =  L  +  NC3  *  4 
DO  211  1=1, NR 
J  =  IR(I) 

211  AA  =  AA  +  MODE(K, J)  *  X1(I+M) 

ENDIF 
T(K)  =  AA 
N  =  N+l 

202  CONTINUE 
WRITE(*,9)A,T(1),  T(2) 

IF  (A.GE.20.0)  GOTO  210 
GOTO  203 

ELSE 

GOTO  201 
ENDIF 

210  CONTINUE 

PRINT*,'  210  ' 

PRINT*, ' 
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FORMAT (3X, F4. 1, 4X, FIS. 6, 4X,  F15. 6) 

END 

FUNCTION  DOT ( NR, A, B) 

DIMENSION  A(l), B(l) 

DOT=0. 

DO  1  1=1, NR 
1  DOT=DOT+A(I)*B(I) 

RETURN 

END 

SUBROUTINE  VADD(N,C1,A,B) 

DIMENSION  A(1),B(1) 

DO  1  1=1, N 

1  A(I)=A(I)+C1*B(I) 

RETURN 

END 

FUNCTION  XNORMl(N, A) 

C  COMPUTES  AN  APPROXIMATION  TO  NORM  OF  A—  NOT  A  BOUND 

DIMENSION  A(N) 

COMMON/MAINA/  NDA.NDAl 

NN=N*NDA 

C1=0. 

TR=A(1) 

IF  (N.EQ.l)  GO  TO  20 
1=2 

DO  10  II=NDA1, NN, NDA 
0=11 

DO  5  00=1, II, NDA 
C1=C1+ABS(A(0)*A(00)) 

5  0=0+1 

TR=TR+A(0) 

10  1=1+1 

TR=TR/F10AT(N) 

DO  15  11=1, NN, N0A1 
15  C1=C1+(A( II )-TR)++2 

20  XN0RM1 =ABS ( TR ) +SQRT ( Cl ) 

RETURN 

END 

FUNCTION  XNORM(N.A) 

C  COMPUTES  AN  APPROXIMATION  TO  NORM  OF  A  —  NOT  A  BOUND 
DIMENSION  A(l) 

COMMON/MAINB/NCOL,  NC0L1 

NN=N*NCOL 

C1=0. 

TR=A(1) 

IF  (N.EQ.l)  GO  TO  20 
1=2 

DO  10  II=NC0L1,NN, NCOL 
0=11 

DO  5  00=1, II, NCOL 
C1=C1+ABS(A(0)*A(00) ) 

5  0=0+1 

TR=TR+A(0) 

10  1=1+1 
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TR=T  R / FLOAT ( N) 

00  15  II=1,NN, NC0L1 
15  Cl=Cl+(A(II)-TR)*-*2 

2«  XNORM= ABS ( TR ) +SQR  T ( Cl ) 

RETURN 

END 

SUBROUTINE  A8GC(A, B,  C, N2, KCC,  NACT,  NSEN, ABG) 
COMMON/MAINB/NCOL, NCOL 1 

mK^«’1)'C(W0L’l)'KCC<NC<x'1,'“e(N^'1) 

CALL  MMUL(B,  KCC, N2, NACT, NACT,  ABG) 

CALL  MMUL(ABG, C, N2, NACT , N2, KCC) 

CALL  UBOAT (A, KCC, N2, N2, ABG) 

RETURN 

END 

SUBROUTINE  MADD(A, B, N,M,  C) 

COMMON/MAINB/NCOL 

R  EAL  A  ( NC  OL ,  NCOL ) ,  B  ( NC  OL ,  NCOL ) ,  C  ( NC  OL ,  NC  OL ) 

INTEGER  N,M 
DO  1  1=1, N 
DO  1  J=1,M 
C(I,0)=A(I, J)+B(I, J) 

1  CONTINUE 

RETURN 
END 

SUBROUTINE  UBOAT(A,B,N,M,C) 

COMMON/MAINB/NCOL 

REAL  A (NCOL, NCOL) , B(NCOL, NCOL), C( NCOL, NCOL) 

INTEGER  N,M 

DO  1  1=1, N 

DO  1  J=1,M 

C(I,  J)=A(I, J)-B(I, J) 

1  CONTINUE 

RETURN 
END 

SUBROUTINE  FINDK(C,M,N,W,G,L,K) 

REAL  UT ( 21 , 21 ) , VP (21,21), QPLUS( 21 , 21 ) 

REAL  SI ( 24) , CPLUS( 21 , 21 ) , QUT (21,21) 

REAL  V(21,21),C(21,21),G(21, 21) 

REAL  K(21, 21), W(48f 48) 

DO  1  1=1, M 
DO  1  J=lfM 

1  UT(I,J)=0.0 
DO  2  1=1, M 

2  UT { I , I ) =1 . 0 
DO  3  1=1, M 
DO  3  J=1,N 

3  VP(I,J)=C(I,J) 

CALL  LSVDF ( VP , 21 , M, N, UT , 21 , M, SI ,  W,  I ER ) 

DO  4  1=1, N 
DO  4  J=1,M 

4  QPLUS(I, J)=0. 0 
DO  5  1=1,  N 


IF(SI(I).GT. 0.M901)  QPLUSfl,  I )=1/SI { I ) 
DO  6  1=1, N 
DO  6  J=1,N 
V(I,J)=VP(I,0) 

CALL  MMUL  ( QPLUS,  UT,  N,  M,  M,  QUT ) 

CALL  MMUL (V, QUT, N,N,M,CPLUS) 

CALL  MMUL(G,CPLUS,L,N,M,K) 

RETURN 

END 


4'D  Matrix  (12x21) 


Mode 

4 

0.004775 

0.004876 

-.004429 

0.000000 

0.004775 

0.004876 

-.004429 

0.000000 

-.013837 

0.003015 

-.002568 

-.013837 

0.003015 

-.002568 

0.009529 

0.009529 

0.000000 

-.009082 

-.013837 

-.009082 

-.013837 

Mode 

5 

0.000001 

0.003619 

0.003618 

-.004642 

0.000001 

-.003618 

-.003619 

0.013451 

-.000002 

0.003619 

0.003618 

-.000002 

-.003618 

-.003619 

0.003620 

-.003617 

0.013452 

-.003620 

-.000002 

0.003617 

-.000002 

Mode 

6 

-.006023 

0.000267 

0.000266 

-.008232 

0.006023 

-.000266 

-.000266 

-.003888 

-.006023 

0.000266 

0.000266 

0.006023 

-.000266 

-.000266 

0.000266 

-.000266 

0.015687 

-.000266 

-.006023 

0.000266 

0.006023 

Mode 

7 

0.000026 

-.001570 

-.001533 

0.002816 

-.000026 

0.001569 

0.001532 

-.006275 

0.000098 

-.001565 

-.001541 

-.000096 

0.001563 

0.001540 

-.001585 

0.001583 

-.006567 

0.001519 

0.000089 

-.001520 

-.000088 

Mode 

12 

-.000452 

-.000859 

-.000903 

0.010104 

0.000432 

0.000869 

0.000927 

-.001977 

-.000531 

-.000880 

-.000909 

0.000536 

0.000891 

0.000928 

-.000822 

0.000824 

-.000105 

0.000970 

-.000570 

-.000940 

0.000577 

Mode 

13 

-.016036 

-.002018 

0.006244 

0.000073 

-.016045 

-.002002 

0.006278 

-.000012 

0.000644 

-.000401 

0.004628 

0.000668 

-.000382 

0.004658 

-.006268 

-.006263 

0.000033 

0.010546 

0.000658 

0.010503 

0.000686 

Mode 

17 

-.000811 

-.000616 

-.000255 

0.000002 

-.000812 

-.000619 

-.000257 

-.000000 

-.000085 

-.000550 

-.000331 

-.000086 

-.000552 

-.000333 

-.000810 

-.000812 

-.000000 

-.000075 

-.000085 

-.000073 

-.000086 

3>"*"d  Matrix  (12x21)  -  (Continued) 


Mode 

21 

-.001659 

-.001165 

0.000520 

-.001659 

0.001649 

0.001152 

-.000527 

0.001115 

0.001645 

-.000847 

0.000291 

-.006164 

0.000840 

-.000291 

-.002043 

0.002030 

-.004348 

-.001337 

0.001618 

0.001328 

-.001622 

Mode 

22 

0.000981 

-.002849 

-.003470 

0.001226 

-.000980 

0.002859 

0.003478 

0.000638 

-.000292 

-.002955 

-.003384 

0.000293 

0.002965 

0.003389 

-.002561 

0.002570 

0.001357 

0.003795 

-.000201 

-.003782 

0.000206 

Mode 

24 

-.019403 

-.005697 

0.005622 

-.017699 

0.019497 

0.005625 

-.005680 

0.001840 

0.002955 

-.003584 

0.004155 

-.002891 

0.003542 

-.004167 

-.011645 

0.011593 

-.008292 

-.011342 

0.002726 

0.011238 

-.002757 

Mode 

28 

-.000000 

-.000004 

0.000001 

-.000003 

-.000003 

-.000005 

0.000001 

-.000000 

-.000008 

0.000001 

-.000004 

-.000008 

-.000000 

-.000003 

0.000004 

0.000003 

0.000000 

-.000007 

-.000009 

-.000009 

-.000009 

Mode 

30 

-.000003 

0.000010 

0.000015 

-.000005 

-.000009 

0.000008 

0.000016 

0.000000 

-.000004 

0.000018 

0.000016 

-.000004 

0.000016 

0.000016 

0.000023 

0.000018 

0.000000 

0.000014 

-.000005 

0.000011 

-.000004 

$TD  Matrix  (33x21) 


Mode  4 

.004775 

.004876 

-.004430 

.004775 

.004876 

-.004430 

-.013837 

.003015 

-.002568 

.003015 

-.002568 

.009529 

-.000000 

-.009082 

-.013837 

-.009082 

Mode  5 

. 000001 

.003619 

.003618 

.000001 

-.003618 

-.003619 

-.000002 

.003619 

.003618 

-.003618 

-.003619 

.003620 

.013452 

-.003620 

-.000002 

.003617 

Mode  6 

-.006023 

.000267 

.000266 

.006023 

-.000266 

-.000266 

-.006023 

.000266 

.000266 

-.000266 

-.000266 

.000267 

.015687 

-.000267 

-.006023 

.000266 

Mode  7 

.000026 

-.001570 

-.001532 

-.000026 

.001569 

.001532 

.000098 

-.001564 

-.001541 

.001563 

.001540 

-.001585 

-.006567 

.001520 

.000090 

-.001520 

Mode  8 

-.000057 

-.002771 

-.002719 

-.000057 

-.002771 

-.002719 

.000049 

-.002762 

-.002730 

-.002762 

-.002730 

-.002799 

-.000000 

-.002693 

.000049 

-.002693 

Mode  9 

.003516 

.000906 

-.000769 

.003515 

.000906 

-.000769 

.000165 

.000572 

-.000434 

.000572 

-.000434 

.001745 

-.000001 

-.001607 

.000165 

-.001607 

Mode  10 

-.005664 

.000499 

.000306 

.005664 

-.000497 

-.000305 

-.006042 

.000463 

.000333 

-.000461 

-.000332 

.000599 

.016056 

-.000213 

-.006038 

.000215 

.000000 

-.000000 

-.013837 

.009529 

-.013837 


-.004642 

.013452 

-.000002 

-.003617 

-.000002 


-.008232 

-.003888 

.006023 

-.000266 

.006023 


.002816 

-.006275 

-.000096 

.001583 

-.000088 


.000002 

-.000000 

.000049 

-.002799 

.000049 


.000000 

.000001 

.000164 

.001745 

.000163 


-.008190 

-.003592 

.006041 

-.000597 

.006038 


$  D  Matrix  (33x21)  - 

(Continued) 

Mode  11 

.001926 

.009437 

.009099 

-.000042 

.001924 

.009438 

.009099 

.000005 

.001279 

.009391 

.009196 

.001274 

.009392 

.009196 

.009638 

.009639 

-.000007 

.001284 

.008962 

.001279 

.008961 

Mode  12 

-.000452 

-.000859 

-.000904 

.010105 

.000432 

.000869 

.000927 

-.001977 

-.000531 

-.000880 

-.000909 

.000536 

.000891 

.000929 

-.000822 

. 000824 

-.000105 

-.000570 

-.000940 

.000577 

.000970 

Mode  13 

-.016036 

-.002018 

.006244 

.000073 

-.016045 

-.002002 

.006278 

-.000012 

. 000644 

-.000401 

.004628 

.000668 

-.000382 

.004658 

-.006268 

-.006263 

.000033 

.000658 

.010503 

.000686 

.010546 

Mode  14 

.000554 

-.001865 

-.001122 

.007250 

-.000518 

.001872 

.001116 

.001287 

.002022 

-.001729 

-.001220 

-.002008 

.001735 

.001220 

-.002262 

.002277 

-.005328 

.002010 

-.000762 

-.002002 

.000748 

Mode  15 

-.006917 

-.001207 

.001357 

-.000168 

-.006931 

-.001238 

.001291 

.000031 

-.001849 

-.000686 

. 000834 

-.001932 

-.000724 

.000773 

-.002469 

-.002481 

-.000105 

-.001852 

.002606 

-.001935 

.002522 

Mode  16 

-.000400 

-.001771 

-.003199 

.014531 

.000377 

.001782 

.003200 

-.002532 

-.003176 

-.002052 

-.003001 

.003124 

.002057 

.  .002995 

-.001045 

.001058 

.007844 

-.003162 

-.003875 

.003123 

.003873 

Mode  17 

-.000911 

-.000616 

-.000255 

.000002 

-.000813 

-.000619 

-.000257 

-.000000 

-.000085 

-.000550 

-.000331 

-.000086 

-.000552 

-.000333 

-.000810 

-.000812 

-.000000 

-.000085 

-.000073 

-.000086 

-.000075 


$  D  Matrix  (33x21)  -  (Continued) 


Mode  18 

.002003 

.000107 

-.000857 

.000002 

.002006 

.000106 

-.000859 

.000000 

.000053 

-.000082 

-.000671 

.000054 

-.000082 

-.000672 

. 000606 

.000606 

.000001 

-.001362 

.000050 

-.001360 

.000051 

Mode  19 

. 000000 

.000000 

-.000000 

-.000001 

.000001 

-.000000 

-.000000 

-.000000 

-.000000 

-.000000 

.000000 

-.000000 

. 000000 
-.000000 

-.000000 

-.000000 

.000000 

Mode  20 

-.000479 

. 000083 

.000025 

-.000018 

-.000475 

. 000088 

.000027 

.000001 

-.000675 

.000089 

.000022 

-.000676 

.000094 

.000025 

.000162 

.000170 

-.000003 

-.000057 

-.000680 

-.000057 

-.000683 

Mode  21 

-.001659 

-.001165 

.000520 

-.001659 

.001650 

.001152 

-.000527 

.001115 

.001645 

-.000847 

.000291 

-.001636 

.000840 

-.000291 

-.002043 

.002030 

-.004348 

-.001337 

.001618 

.001328 

-.001621 

Mode  22 

.000981 

-.002849 

-.003470 

.001226 

-.000980 

.002859 

.003478 

.000638 

-.000292 

-.002955 

-.003384 

.000293 

.002965 

.003389 

-.002561 

.002570 

.001357 

-.000201 

-.003782 

.000206 

.003795 

Mode  23 

-.000486 

-.021293 

-.020538 

-.005602 

.000464 

.021322 

.020568 

. 003478 

.000061 

-.021077 

-.020851 

-.000042 

.021124 

.020881 

-.021704 

.021719 

.001631 

.020813 

.000545 

-.020744 

-.000513 

Mode  24 

-.019403 

-.005697 

.005622 

-.017699 

$  D  Matrix  (33x21)  - 

(Continued) 

Mode  25 

.000001 

.000000 

-.000000 

.000001 

-.000001 

-.000000 

.000000 

-.000000 

-.000000 

.000000 

-.000000 

.000000 

-.000000 

.000000 

.000000 

-.000000 

.000000 

-.000000 

-.000000 

.000000 

.000000 

Mode  26 

-.000001 

-.000001 

-.000001 

-.000001 

-.000004 

-.000002 

-.000001 

-.000000 

-.000012 

.000001 

-.000003 

-.000011 

.000000 

-.000003 

. 000006 

.000003 

. 000000 

-.000012 

-.000008 

-.000012 

-.000007 

Mode  27 

-.000000 

. 000000 

-.000000 

. 000000 

.000000 

. 000000 

-.000000 

.000000 

. 000000 

-.000000 

.000000 

. 000000 

.000000 

.000000 

-.000000 

-.000000 

-.000000 

. 000000 

.000000 

. 000000 

.000000 

Mode  28 

-.000000 

-.000004 

.000001 

-.000003 

-.000003 

-.000005 

.000001 

-.000000 

-.000008 

.000001 

-.000004 

-.000008 

-.000000 

-.000003 

.000004 

.000003 

.000000 

-.000009 

-.000009 

-.000009 

-.000007 

Mode  29 

-.000043 

-.000006 

. 000008 

.  000000 

-.000042 

-.000006 

.  000008 

-.000000 

-.000004 

-.000008 

.000006 

- . 000004 

-.000008 

.000006 

-.000027 

-.000027 

-.000000 

-.000003 

. 000024 

-.000003 

. 000024 

Mode  30 

-.000003 

.000010 

.000015 

-.000005 

-.000009 

.000008 

.000016 

.000000 

-.000004 

.000018 

.000016 

-.000004 

.000016 

.000016 

.000023 

.000018 

.000000 

-.000005 

.000011 

-.000004 

.000014 

Mode  31 

-.000195 

-.000069 

.000030 

-.000195 

.000197 

.000069 

-.000030 

.000001 

.000002 

-.000051 

.000019 

-.000001 

.000052 

-.000019 

-.000124 

.000125 

-.000006 

.000001 

.000085 

-.000001 

-.000086 


Matrix  (33x21)  -  (Continued) 


Mode 

32 

-.000070 

.000105 

.000112 

-.000011 

.000071 

-.000105 

-.000111 

.000001 

. 000004 

.000107 

.000137 

-.000002 

-.000107 

-.000136 

.000081 

-.000080 

-.000020 

-.000170 

.000005 

.000169 

-.000006 

Mode 

33 

-.001316 

-.002646 

.001407 

-.000041 

-.001433 

-.003269 

.001633 

.000746 

.001072 

.008585 

.006652 

.002688 

.008136 

.006664 

.017513 

.016781 

.000026 

-.000294 

.000412 

-.000583 

.001442 

Mode 

34 

-.011539 

-.001697 

.003121 

.002881 

.011322 

-.000778 

-.003301 

.000771 

-.000215 

-.011488 

-.007493 

.003215 

.010932 

.006959 

-.021277 

.023013 

.001047 

-.014854 

.001293 

.010261 

-.003702 

Mode 

35 

-.005960 

.007119 

.007164 

.003841 

.000284 

.009320 

.007115 

-.002497 

-.002435 

-.000120 

.002235 

-.006592 

.004583 

.001767 

-.008845 

.000859 

-.000228 

.002261 

.018397 

-.000966 

.012288 

Mode 

36 

-.010398 

-.003798 

-.005374 

.014434 

.010019 

-.000317 

.001925 

-.000109 

-.000137 

-.004307 

.003640 

.001805 

.003038 

-.004431 

-.010575 

.011795 

-.001445 

-.010171 

-.000139 

.004044 

.000032 

t 


I 


Bib! iography 


I* 


h 


El  * 


fcs 


N 


L< 


M 


1.  Aldridge,  2nd  Lt  Edward  S.  Decentralized  Control  of  a  Large  Space 
Structure  as  Applied  to  the  CSDL  2  Model.  MS  Thesis,  GA/AA/82D-1. 
School  of  Engineering,  Air  Force  Institute  of  Technology  (AU),  Wright- 
Patterson  AFB  OH,  December  1982  (AD-A124702) . 

2.  Calico,  Robert  A.  Jr.  and  W.  T.  Miller.  "Decentralized  Control  of  a 
Flexible  Spacecraft,"  AIAA/AAS  Astrodynamics  Conference.  AIAA-82-1414. 
San  Diego,  CA  August  9-ll,  1982. 

3.  D'Azzo,  John  J.  and  Constantine  H.  Houpis.  Linear  Control  System 
Analysis  and  Design:  Conventional  and  Modern.  (Second  Edition) . 

New  York:  McGraw-Hill  Inc.,  1981. 

4.  Henderson,  Timothy.  VCQSS  Design  Model  arid  Disturbances.  Revised 
Report  for  CSDL  Publications  C-534/,  Charles  Stark  Draper  Laboratory, 
Inc.,  Cambridge  Massachusetts,  November  12,  1981. 

5.  - .  Active  Control  of  Space  Structures  (ACOSS)  Model  2.  Work 

performed  by  the  Charles  Stark  Draper  Laboratory,  Inc.  (CSDL)  under 
Contract  F30602-80-C-0096  as  part  of  the  ACOSS  program  for  the  Rome 
Air  Development  Center,  September  1981. 

6.  Hess,  Jeffrey  W.  Unpublished  notes  on  Decentralized  Control  Using  Di¬ 
rect  Output  Feedback.  School  of  Engineering,  Air  Force  Institute  of 
Technology,  Wright-Patterson  AFB  OH,  1982. 

7.  Kleinman,  D.  L.  A  Description  of  Computer  Programs  for  Use  in  Linear 
Systems  Studies.  The  University  of  Connecticut  School  of  Engineering, 
TR-77-2,  Storrs,  Connecticut,  July  1977. 

8.  Laub,  Alan  J.  "A  Schur  Method  for  Solving  Algebraic  Riccati  Equa¬ 
tions,"  IEEE  Transactions  on  Automatic  Control,  Vol .  AC-24,  No.  6, 
p.  913. 

9.  Maybeck,  Peter  S.  Stochastic  Models,  Estimation,  and  Control ,  Vols. 
1-3.  New  York:  Academic  Press,  1979,  1982. 

10.  Meirovitch,  L.  et  al_.  "A  Comparison  of  Control  Techniques  for  Large 
Flexible  Systems,"  Journal  of  Guidance,  Control  and  Dynamics,  Vol.  6, 
No.  4,  July-August  1983. 

11.  Meirovitch,  L.  and  H.  Oz.  "Computational  Aspects  of  the  Control  of 
Large  Flexible  Structures,"  Proceedings  of  18th  IEEE  Conference  on 
Decision  and  Control .  December  12-14, 1979. 

12.  - .  "Modal -Space  Control  of  Distributed  Gyroscopic  Systems," 

Journal  of  Guidance  and  Control ,  Vol.  3,  No.  2,  p.  140,  1978. 


125 


13.  Miller,  William  T.  ''Decentral  ized  Control  of  Large  Space  Structures. " 
MS  Thesis,  School  of  Engineering,  Air  Force  Institute  of  Technology, 
Wright-Patterson  AFB  OH,  December  1981. 

14.  RADC-TR-83-56.  Modifications  to  ACOSS  Model  #2  Design.  A  Study  con¬ 
ducted  by  Charles  Stark  Draper  Lab,  Inc.  for  Rome  Air  Development 
Center,  March  1983. 

15.  RADC-TR-82-21.  ACOSS  FIVE  (Active  Control  of  Space  Structures)  Phase 
IA.  Study  conducted  by  Lockheed  Missiles  and  Space  Company,  Inc.  for 
the  Rome  Air  Development  Center,  March  1982  (AD-A116655) . 

16.  RADC-TR-81-242.  ACOSS  EIGHT  (Active  Control  of  Space  Structures) 

Phase  II.  Study  conducted  by  TRW  for  Rome  Air  Development  Center, 
September  1981  (AD-A108629) . 

17.  RADC-TR-80-78.  ACOSS  FOUR  (Active  Control  of  Space  Structures)  Theory, 
Vol  I  (of  two).  Study  conducted  by  Charles  Stark  Draper  Lab,  Inc. 

for  Rome  Air  Development  Center,  April  1980  (AD-A085654) . 

18.  RADC-TR-79-268.  Actively  Controlled  Structures  Theory.  Theory  of 
Design  Methods,  Vol  I  (of  two).  Study  conducted  by  Charles  Stark 
Draper  Laboratory,  Inc.  for  Rome  Air  Development  Center,  November 
1979  (AD-8044535) . 

19.  Salzwedel,  Horst  and  J.  H.  Vincent.  "Modeling,  Identification,  and 
Control  of  Flexible  Aircraft,"  AIAA  Atmospheric  Flight  Mechanics  Con¬ 
ference.  Gatl inburg,  TN,  August  15-17,  1983. 

20.  Sandell,  Nils  et  al_.  "Survey  of  Decentralized  Control  Methods  for 
Large  Scale  Systems,"  IEEE  Transactions  on  Automatic  Control,  Vol. 
AC-23,  No.  2,  p.  108. 

21.  Sesak,  J.  R.  et^  "Flexible  Spacecraft  Control  by  Model  Error  Sen¬ 
sitivity  Suppression,"  The  Journal  of  the  Astronautical  Sciences, 

Vol.  27,  No.  2,  pp.  131-156,  April -June  1974"! 

22.  Strang,  Gilbert.  Linear  Algebra  and  Its  Application.  New  York:  Aca¬ 
demic  Press,  1976. 


David  Vernon  Thyfault  is  a  native  of  Colorado.  He  received  his  Bach¬ 


elor  of  Science  in  Aerospace  Engineering  Sciences  from  the  University  of 
Colorado  in  May  1974  and  was  contnissioned  as  a  second  lieutenant  in  the 
United  States  Air  Force  through  the  Officer  Training  School,  Lackland  AFB, 
Texas  in  September  1974.  He  served  as  a  space  vehicle  engineer  with  the 
6595th  Space  Test  Group,  Special  Space  Projects  Branch,  Vandenberg  AFB, 
California.  While  at  Vandenberg,  he  participated  in  design,  test,  and 
launch  activities  for  the  Defense  Meteorological  Satellite  Program,  Space 
Test  Program  vehicles  P72-2,  P76-5,  P80-1,  and  P80-2,  and  the  NASA  meteor¬ 
ological  series  TIROS/NOAA  program.  His  following  assignment  was  as  a 
space  laser  project  officer  for  the  Air  Force  Weapons  Laboratory,  Advanced 
Beam  Control  Branch,  Kirtland  AFB,  New  Mexico,  where  he  evaluated  control 
system  concepts  for  high  energy  laser  programs. 


Permanent  address:  1849  Mariposa  Avenue 


Boulder,  Colorado  80302 


UNCLASSIFIED 


security  classification  of  this  page 


REPORT  DOCUMENTATION  PAGE 


lb.  RESTRICTIVE  MARKINGS 


2a.  SECURITY  CLASSIFICATION  AUTHORITY 


2b.  DECLASSIFICATION/DOWNGRAOING  SCHEDULE 


4.  PERFORMING  ORGANIZATION  REPORT  NUMBERIS) 

AFIT/GA/AA/33-D-S 


3.  DISTRIBUTION/AVAILABILITY  OF  REPORT 

Approved  for  public  release; 
distribution  unlimited. 


5.  MONITORING  ORGANIZATION  REPORT  NUMBERIS) 


6a  NAME  OF  PERFORMING  ORGANIZATION  |6b.  OFFICE  SYMBOL  7a.  NAME  OF  MONITORING  ORGANIZATION 

I  (If  applicable! 

School  of  Engineering  | 


6c.  AOORESS  (City.  State  and  ZIP  Code! 


7b.  ADDRESS  (City,  State  and  ZIP  Code I 


Air  Force  Institute  of  Technology 
Wripht-Patterson  AFB,  Ohio  45433 


8a.  NAME  OF  FUNOING/SPONSORiNG  8b.  OFFICE  SYMBOL  9.  PROCUREMENT  INSTRUMENT  IDENTIFICATION  NUMBER 

organization  Ai r  Force  Flight  a r applicable) 

Dynamics  Laboratory  "  FIGCA 


8c.  ADDRESS  (City.  State  and  ZIP  Code) 

Wriaht-Patterson  AFB,  Ohio  45433 


IV  TITLE  ( Include  Security  Classification) 

See  Box  19 


'J.  PERSONAL  AUTHOR(S) 

David  V.  Thyfault,  Capt. ,  USAF 


13a  TYPE  OF  REPORT  13b.  TIME  COVERED 

MS  Thesis  prom _ 


16.  supplementary  notation 


1 10.  SOURCE  OP  FUNDING  NOS. 

PROGRAM 

project 

TASK 

WORK  UNIT 

element  no. 

NO. 

NO 

NO 

14.  DATE  OF  REPORT  (Yr  .  Mo..  Day)  15.  PAGE  COUNT 

1983  December  .  .  mi  137 


BS 


3  vf.') 


p,  f,  .  .  «  P*  r. '  Tr "clcpmeaf 

_ _  _  _ HiT  fitrrn  I _ '1— ^  C  f  ‘J  I  _ 


17.  COSATI  CODES _  18.  SUBJECT  TERMS  <CorMH&'*frti*6krk  fhtocedory~aKd  identify  by  block  number ) 

GROUP _ SUB.  GR.  ^JClosed  Loop  System,  Decentralization,  Equations  of  State, 

02  IControl  Theory,  Actuators,  Riccati  Eauation,  Feedback, 

Decoupling  (Interaction),  Stability,  Control  Systems 


18.  ABSTRACT  i Continue  on  reverse  if  necessary  and  identify  by  block  number ; 

Title:  DECENTRALIZED  CONTROL  OF  A  LARGE  SPACE 
STRUCTURE  USING  DIRECT  OUTPUT  FEEDBACK 

Thesis  Chairman:  Dr.  Robert  A.  Calico,  Jr. 


•  •a.  oistribution/availability  of  abstract 
UNCLASSIFIED/UNLIMITED  QE  SAME  AS  RPT.  □  OTIC  USERS  □ 


21  ABSTRACT  SECURITY  CLASSIFICATION 

UNCLASSIFIED 


22a.  NAME  OF  RESPONSIBLE  INDIVIDUAL 

22b.  TELEPHONE  NUMBER 
(Include  Area  Code) 

22c  OFFICE  symbol 

Robert  A.  Calico,  Jr.,  Professor 

(513)  255-2104 

AFIT/ENY 

DD  FORM  1473, 83  APR 


EOITION  OF  1  JAN  73  IS  OBSOLETE. 


UNCLASSIFIED _ 

security  classification  of  this  page 


UNCLASSIFIED 


security  classification  of  this  page 


Direct  output  feedback  control  methods  are  used  to  develop  a 
multiple-input  multiple-output  controller.  The  controller  is  applied 
to  the  ACOSS  2  Model  (CSDL  2).  NASTRAN  is  employed  to  generate  modal 
approximations  of  the  model.  The  first  36  modes  are  utilized  and  im¬ 
plemented  in  the  controller.  The  control  problem  is  formulated  in 
state  space  form  and  direct  output  feedback  is  implemented.  The  state 
is  represented  as  modal  amplitudes  and  rates.  System  outputs  are 
obtained  by  rate  sensors  and  control  is  applied  by  point  force  actuators. 
Since  the  output  matrices  are  not  of  full  rank,  their  generalized 
inverses  are  obtained.  These  and  steady  full -state  feedback  optimal 
gains  are  used  to  determine  suboptimal  feedback  gains.  Control  and 
observation  spillover  are  eliminated  using  singular  value  decomposition. 
Decentralized  control  is  accomplished  using  three  and  four  controllers. 
Conditions  for  which  the  stability  of  the  model  is  assured  are  developed. 
Sensitivity  to  uncontrolled  residual  modes  was  examined.  Direct  output 
feedback  response  was  compared  to  that  obtained  using  optimal  ful 1  - 
state  feedback.  Full  controller  decoupling  was  achieved  and  stability 
maintained. 


UNCLASSIFIED 


SECURITY  1*1  ASSlFlCATION  of  this  cage 


